From ed091655db0a5411ce6dbcc98bdfac07f1431169 Mon Sep 17 00:00:00 2001
From: "AMATERASU\\geier" <geier@irmb.tu-bs.de>
Date: Tue, 15 Jun 2021 14:23:52 +0200
Subject: [PATCH] Two phase fields implemented but no energy equation.
---
apps/cpu/Multiphase/Multiphase.cpp | 2 +-
apps/cpu/Multiphase/MultiphaseGeier.cfg | 60 ++
.../WriteMultiphaseQuantitiesCoProcessor.cpp | 52 +-
.../MultiphaseScratchCumulantLBMKernel.cpp | 522 ++++++++--
...tiphaseTwoPhaseFieldsCumulantLBMKernel.cpp | 940 +++++++++++++++---
...ultiphaseTwoPhaseFieldsCumulantLBMKernel.h | 5 +
...ultiphaseInitDistributionsBlockVisitor.cpp | 31 +-
7 files changed, 1397 insertions(+), 215 deletions(-)
create mode 100644 apps/cpu/Multiphase/MultiphaseGeier.cfg
diff --git a/apps/cpu/Multiphase/Multiphase.cpp b/apps/cpu/Multiphase/Multiphase.cpp
index deb2845f4..a22777658 100644
--- a/apps/cpu/Multiphase/Multiphase.cpp
+++ b/apps/cpu/Multiphase/Multiphase.cpp
@@ -137,7 +137,7 @@ void run(string configname)
fctF2.SetExpr("vy1");
fctF2.DefineConst("vy1", uLB);
- double startTime = 500;
+ double startTime = 30;
SPtr<BCAdapter> velBCAdapterF1(new MultiphaseVelocityBCAdapter(true, false, false, fctF1, phiH, 0.0, startTime));
SPtr<BCAdapter> velBCAdapterF2(new MultiphaseVelocityBCAdapter(true, false, false, fctF2, phiH, startTime, endTime));
diff --git a/apps/cpu/Multiphase/MultiphaseGeier.cfg b/apps/cpu/Multiphase/MultiphaseGeier.cfg
new file mode 100644
index 000000000..bcabb0684
--- /dev/null
+++ b/apps/cpu/Multiphase/MultiphaseGeier.cfg
@@ -0,0 +1,60 @@
+#pathname = E:/Multiphase/HesamCodeWithCumulantsDensRatio
+#pathname = E:/Multiphase/HesamCodeWithCumulantsQuartic
+#pathname = E:/Multiphase/HesamCode
+pathname = E:/Multiphase/HesamCodeCumulantTubeFilter
+pathGeo = C:/Users/geier/Documents/VirtualFluids_dev_Kostya/apps/cpu/Multiphase/backup
+geoFile=tubeTransformed.stl
+#geoFile = JetBreakup2.ASCII.stl
+numOfThreads = 4
+availMem = 10e9
+
+#Grid
+
+#boundingBox = -1.0 121.0 0.5 629.0 -1.0 121.0 #(Jet Breakup) (Original with inlet length)
+#boundingBox = -60.5 60.5 -1.0 -201.0 -60.5 60.5 #(Jet Breakup2) (Original without inlet length)
+#blocknx = 22 20 22
+
+#boundingBox = -60.5 60.5 -1.0 -21.0 -60.5 60.5 #(Jet Breakup2) (Original without inlet length)
+#boundingBox = -60.5 60.5 -21.0 -1.0 -60.5 60.5 #(Jet Breakup2) (Original without inlet length)
+#blocknx = 22 20 22
+
+
+#dx = 0.5
+
+#boundingBox = 6.0e-3 46.0e-3 -5e-3 5e-3 -5e-3 5e-3
+boundingBox = 6.0e-3 16.0e-3 -5e-3 5e-3 -5e-3 5e-3
+blocknx = 60 60 60 #20 20 20
+
+dx = 1.66666666667e-4
+
+refineLevel = 0
+
+#Simulation
+uLB =0.005# 0.0000005 #inlet velocity
+uF2 = 0.0001
+Re = 10
+nuL =1e-6#1e-2# 1.0e-5 #!1e-2
+nuG =1e-6#1e-2# 1.16e-4 #!1e-2
+densityRatio = 10000#1000#1000 #30
+sigma =0# 1e-4 #4.66e-3 #surface tension 1e-4 ./. 1e-5
+interfaceThickness = 5
+radius = 615.0 (Jet Breakup)
+contactAngle = 110.0
+gravity = 0.0
+#gravity = -5.04e-6
+phi_L = 0.0
+phi_H = 1.0
+Phase-field Relaxation = 0.6
+Mobility = 0.1 #0.02 # 0.01 ./. 0.08, fine correction of Phase-field Relaxation parameter, to activate it need to change in kernel tauH to tauH1
+
+
+logToFile = false
+
+newStart = true
+restartStep = 100000
+
+cpStart = 100000
+cpStep = 100000
+
+outTime = 100
+endTime = 200000000
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
index 73034d889..70d2f2b6c 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
@@ -155,6 +155,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
datanames.push_back("Vy");
datanames.push_back("Vz");
datanames.push_back("P1");
+ datanames.push_back("Phi2");
data.resize(datanames.size());
@@ -162,10 +163,12 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
SPtr<DistributionArray3D> distributionsF = kernel->getDataSet()->getFdistributions();
SPtr<DistributionArray3D> distributionsH = kernel->getDataSet()->getHdistributions();
+ SPtr<DistributionArray3D> distributionsH2 = kernel->getDataSet()->getH2distributions();
SPtr<PhaseFieldArray3D> divU = kernel->getDataSet()->getPhaseField();
LBMReal f[D3Q27System::ENDF + 1];
LBMReal phi[D3Q27System::ENDF + 1];
+ LBMReal phi2[D3Q27System::ENDF + 1];
LBMReal vx1, vx2, vx3, rho, p1, beta, kappa;
LBMReal densityRatio = kernel->getDensityRatio();
@@ -202,6 +205,8 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
CbArray3D<int> nodeNumbers((int)maxX1, (int)maxX2, (int)maxX3, -1);
CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField(
new CbArray3D<LBMReal, IndexerX3X2X1>(maxX1, maxX2, maxX3, -999.0));
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField2(
+ new CbArray3D<LBMReal, IndexerX3X2X1>(maxX1, maxX2, maxX3, -999.0));
for (int ix3 = minX3; ix3 < maxX3; ix3++) {
for (int ix2 = minX2; ix2 < maxX2; ix2++) {
@@ -211,8 +216,17 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
(*phaseField)(ix1, ix2, ix3) =
((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
(((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
- ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
- ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+ ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
+ ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+ if (distributionsH2) {
+ distributionsH2->getDistribution(f, ix1, ix2, ix3);
+ (*phaseField2)(ix1, ix2, ix3) =
+ ((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
+ (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
+ ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
+ ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+ }
+ else { (*phaseField2)(ix1, ix2, ix3) = 999.0; }
}
}
}
@@ -244,6 +258,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
float(worldCoordinates[2])));
phi[REST] = (*phaseField)(ix1, ix2, ix3);
+ phi2[REST] = (*phaseField2)(ix1, ix2, ix3);
if ((ix1 == 0) || (ix2 == 0) || (ix3 == 0)) {
dX1_phi = 0.0;
@@ -284,6 +299,38 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
dX2_phi = 0.0 * gradX2_phi(phi);
dX3_phi = 0.0 * gradX3_phi(phi);
mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi(phi);
+
+ //phi2[E] = (*phaseField2)(ix1 + DX1[E], ix2 + DX2[E], ix3 + DX3[E]);
+ //phi2[N] = (*phaseField2)(ix1 + DX1[N], ix2 + DX2[N], ix3 + DX3[N]);
+ //phi2[T] = (*phaseField2)(ix1 + DX1[T], ix2 + DX2[T], ix3 + DX3[T]);
+ //phi2[W] = (*phaseField2)(ix1 + DX1[W], ix2 + DX2[W], ix3 + DX3[W]);
+ //phi2[S] = (*phaseField2)(ix1 + DX1[S], ix2 + DX2[S], ix3 + DX3[S]);
+ //phi2[B] = (*phaseField2)(ix1 + DX1[B], ix2 + DX2[B], ix3 + DX3[B]);
+ //phi2[NE] = (*phaseField2)(ix1 + DX1[NE], ix2 + DX2[NE], ix3 + DX3[NE]);
+ //phi2[NW] = (*phaseField2)(ix1 + DX1[NW], ix2 + DX2[NW], ix3 + DX3[NW]);
+ //phi2[TE] = (*phaseField2)(ix1 + DX1[TE], ix2 + DX2[TE], ix3 + DX3[TE]);
+ //phi2[TW] = (*phaseField2)(ix1 + DX1[TW], ix2 + DX2[TW], ix3 + DX3[TW]);
+ //phi2[TN] = (*phaseField2)(ix1 + DX1[TN], ix2 + DX2[TN], ix3 + DX3[TN]);
+ //phi2[TS] = (*phaseField2)(ix1 + DX1[TS], ix2 + DX2[TS], ix3 + DX3[TS]);
+ //phi2[SW] = (*phaseField2)(ix1 + DX1[SW], ix2 + DX2[SW], ix3 + DX3[SW]);
+ //phi2[SE] = (*phaseField2)(ix1 + DX1[SE], ix2 + DX2[SE], ix3 + DX3[SE]);
+ //phi2[BW] = (*phaseField2)(ix1 + DX1[BW], ix2 + DX2[BW], ix3 + DX3[BW]);
+ //phi2[BE] = (*phaseField2)(ix1 + DX1[BE], ix2 + DX2[BE], ix3 + DX3[BE]);
+ //phi2[BS] = (*phaseField2)(ix1 + DX1[BS], ix2 + DX2[BS], ix3 + DX3[BS]);
+ //phi2[BN] = (*phaseField2)(ix1 + DX1[BN], ix2 + DX2[BN], ix3 + DX3[BN]);
+ //phi2[BSW] = (*phaseField2)(ix1 + DX1[BSW], ix2 + DX2[BSW], ix3 + DX3[BSW]);
+ //phi2[BSE] = (*phaseField2)(ix1 + DX1[BSE], ix2 + DX2[BSE], ix3 + DX3[BSE]);
+ //phi2[BNW] = (*phaseField2)(ix1 + DX1[BNW], ix2 + DX2[BNW], ix3 + DX3[BNW]);
+ //phi2[BNE] = (*phaseField2)(ix1 + DX1[BNE], ix2 + DX2[BNE], ix3 + DX3[BNE]);
+ //phi2[TNE] = (*phaseField2)(ix1 + DX1[TNE], ix2 + DX2[TNE], ix3 + DX3[TNE]);
+ //phi2[TNW] = (*phaseField2)(ix1 + DX1[TNW], ix2 + DX2[TNW], ix3 + DX3[TNW]);
+ //phi2[TSE] = (*phaseField2)(ix1 + DX1[TSE], ix2 + DX2[TSE], ix3 + DX3[TSE]);
+ //phi2[TSW] = (*phaseField2)(ix1 + DX1[TSW], ix2 + DX2[TSW], ix3 + DX3[TSW]);
+
+ // mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi(phi);
+
+
+
}
distributionsF->getDistribution(f, ix1, ix2, ix3);
@@ -361,6 +408,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
data[index++].push_back(vx2);
data[index++].push_back(vx3);
data[index++].push_back(p1);
+ data[index++].push_back(phi2[REST]);
}
}
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
index 505007e0c..77695179c 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
@@ -39,6 +39,7 @@
#include "DataSet3D.h"
#include "LBMKernel.h"
#include <cmath>
+#include <iostream>
#define PROOF_CORRECTNESS
@@ -163,6 +164,16 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, 0.0));
+ /////For velocity filter
+
+ //CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr velocityX(
+ // new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, 0.0));
+ //CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr velocityY(
+ // new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, 0.0));
+ //CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr velocityZ(
+ // new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, 0.0));
+
+
for (int x3 = 0; x3 <= maxX3; x3++) {
for (int x2 = 0; x2 <= maxX2; x2++) {
for (int x1 = 0; x1 <= maxX1; x1++) {
@@ -207,7 +218,63 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
// (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) +
// (mfbaa + mfbac + mfbca + mfbcc) + (mfabb + mfcbb) +
// (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
+
+ ///Velocity filter
+
+
+ LBMReal rhoH = 1.0;
+ LBMReal rhoL = 1.0 / densityRatio;
+
+ LBMReal rhoToPhi = (rhoH - rhoL) / (phiH - phiL);
+
+
+ LBMReal rho = rhoH + rhoToPhi * ((*phaseField)(x1, x2, x3) - phiH);
+
+ mfbbc = (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3) / rho * c3;
+ mfbcb = (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) / rho * c3;
+ mfccb = (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3) / rho * c3;
+ mfacb = (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3) / rho * c3;
+ mfcbb = (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) / rho * c3;
+ mfcbc = (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3) / rho * c3;
+ mfabc = (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3) / rho * c3;
+ mfbcc = (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3) / rho * c3;
+ mfbac = (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3) / rho * c3;
+ mfccc = (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3) / rho * c3;
+ mfacc = (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3) / rho * c3;
+ mfcac = (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3) / rho * c3;
+ mfaac = (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3) / rho * c3;
+
+ mfabb = (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3) / rho * c3;
+ mfbab = (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3) / rho * c3;
+ mfbba = (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p) / rho * c3;
+ mfaab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3) / rho * c3;
+ mfcab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3) / rho * c3;
+ mfaba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p) / rho * c3;
+ mfcba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p) / rho * c3;
+ mfbaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p) / rho * c3;
+ mfbca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p) / rho * c3;
+ mfaaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p) / rho * c3;
+ mfcaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p) / rho * c3;
+ mfaca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) / rho * c3;
+ mfcca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) / rho * c3;
+
+ mfbbb = (*this->zeroDistributionsF)(x1, x2, x3) / rho * c3;
+
+ //(*velocityX)(x1, x2, x3) = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+ // (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
+ // (mfcbb - mfabb)) ;
+ //(*velocityY)(x1, x2, x3) = ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+ // (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
+ // (mfbcb - mfbab)) ;
+ //(*velocityZ)(x1, x2, x3) = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+ // (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
+ // (mfbbc - mfbba)) ;
+
+
+
+
}
+ else { (*phaseField)(x1, x2, x3) = 0; }
}
}
}
@@ -218,6 +285,10 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
for (int x3 = minX3; x3 < maxX3; x3++) {
for (int x2 = minX2; x2 < maxX2; x2++) {
for (int x1 = minX1; x1 < maxX1; x1++) {
+
+ //for (int x3 = minX3+1; x3 < maxX3-1; x3++) {
+ // for (int x2 = minX2+1; x2 < maxX2-1; x2++) {
+ // for (int x1 = minX1+1; x1 < maxX1-1; x1++) {
if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3)) {
int x1p = x1 + 1;
int x2p = x2 + 1;
@@ -245,7 +316,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
//-1 0 1
findNeighbors(phaseField, x1, x2, x3);
-
+ //// reading distributions here appears to be unnecessary!
LBMReal mfcbb = (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3);
LBMReal mfbcb = (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3);
LBMReal mfbbc = (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3);
@@ -284,6 +355,62 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
LBMReal dX2_phi = gradX2_phi();
LBMReal dX3_phi = gradX3_phi();
+ //LBMReal dX1_phi = 3.0*((
+ // WEIGTH[TNE]*((((*phaseField)(x1 + 1, x2+1, x3+1)- (*phaseField)(x1 - 1, x2 - 1, x3 - 1))+ ((*phaseField)(x1 + 1, x2 - 1, x3 + 1) - (*phaseField)(x1 - 1, x2 + 1, x3 - 1)))
+ // + (((*phaseField)(x1 + 1, x2 - 1, x3 - 1) - (*phaseField)(x1 - 1, x2 + 1, x3 + 1)) + ((*phaseField)(x1 + 1, x2 + 1, x3 - 1) - (*phaseField)(x1 - 1, x2 - 1, x3 + 1))))
+ // +WEIGTH[NE]* ((((*phaseField)(x1 + 1, x2 + 1, x3) - (*phaseField)(x1 - 1, x2 - 1, x3)) + ((*phaseField)(x1 + 1, x2 - 1, x3) - (*phaseField)(x1 - 1, x2 + 1, x3 )))
+ // + (((*phaseField)(x1 + 1, x2, x3 - 1) - (*phaseField)(x1 - 1, x2, x3 + 1)) + ((*phaseField)(x1 + 1, x2, x3 + 1) - (*phaseField)(x1 - 1, x2, x3 - 1)))))
+ // +WEIGTH[N]*((*phaseField)(x1 + 1, x2, x3 ) - (*phaseField)(x1 - 1, x2, x3))
+ // );
+ ////if (dX1_phi != NdX1_phi) {std::cout<<dX1_phi<<" "<< NdX1_phi<<std::endl;}
+
+ //LBMReal dX2_phi = 3.0 * ((
+ // WEIGTH[TNE] * ((((*phaseField)(x1 + 1, x2 + 1, x3 + 1) - (*phaseField)(x1 - 1, x2 - 1, x3 - 1)) + ((*phaseField)(x1 -1, x2 + 1, x3 + 1) - (*phaseField)(x1 + 1, x2 - 1, x3 - 1)))
+ // + (((*phaseField)(x1 - 1, x2 + 1, x3 - 1) - (*phaseField)(x1 + 1, x2 - 1, x3 + 1)) + ((*phaseField)(x1 + 1, x2 + 1, x3 - 1) - (*phaseField)(x1 - 1, x2 - 1, x3 + 1))))
+ // + WEIGTH[NE] * ((((*phaseField)(x1 + 1, x2 + 1, x3) - (*phaseField)(x1 - 1, x2 - 1, x3)) + ((*phaseField)(x1 - 1, x2 + 1, x3) - (*phaseField)(x1 + 1, x2 - 1, x3)))
+ // + (((*phaseField)(x1, x2+1, x3 - 1) - (*phaseField)(x1 , x2-1, x3 + 1)) + ((*phaseField)(x1 , x2+1, x3 + 1) - (*phaseField)(x1 , x2-1, x3 - 1)))))
+ // + WEIGTH[N] * ((*phaseField)(x1 , x2+1, x3) - (*phaseField)(x1 , x2-1, x3))
+ // );
+
+ //LBMReal dX3_phi = 3.0 * ((
+ // WEIGTH[TNE] * ((((*phaseField)(x1 + 1, x2 + 1, x3 + 1) - (*phaseField)(x1 - 1, x2 - 1, x3 - 1)) + ((*phaseField)(x1 - 1, x2 + 1, x3 + 1) - (*phaseField)(x1 + 1, x2 - 1, x3 - 1)))
+ // + (((*phaseField)(x1 - 1, x2 - 1, x3 + 1) - (*phaseField)(x1 + 1, x2 + 1, x3 - 1)) + ((*phaseField)(x1 + 1, x2 - 1, x3 + 1) - (*phaseField)(x1 - 1, x2 + 1, x3 - 1))))
+ // + WEIGTH[NE] * ((((*phaseField)(x1 + 1, x2, x3+1) - (*phaseField)(x1 - 1, x2, x3-1)) + ((*phaseField)(x1 - 1, x2, x3+1) - (*phaseField)(x1 + 1, x2, x3-1)))
+ // + (((*phaseField)(x1, x2 - 1, x3 + 1) - (*phaseField)(x1, x2 + 1, x3 - 1)) + ((*phaseField)(x1, x2 + 1, x3 + 1) - (*phaseField)(x1, x2 - 1, x3 - 1)))))
+ // + WEIGTH[N] * ((*phaseField)(x1, x2, x3+1) - (*phaseField)(x1, x2, x3-1))
+ // );
+
+ ///////////////////////////////////////
+
+ //LBMReal dX1_phi2 = 1.5 * ((
+ // WEIGTH[TNE] * ((((*phaseField)(x1 + 2, x2 + 2, x3 + 2) - (*phaseField)(x1 - 2, x2 - 2, x3 - 2)) + ((*phaseField)(x1 + 2, x2 - 2, x3 + 2) - (*phaseField)(x1 - 2, x2 + 2, x3 - 2)))
+ // + (((*phaseField)(x1 + 2, x2 - 2, x3 - 2) - (*phaseField)(x1 - 2, x2 + 2, x3 + 2)) + ((*phaseField)(x1 + 2, x2 + 2, x3 - 2) - (*phaseField)(x1 - 2, x2 - 2, x3 + 2))))
+ // + WEIGTH[NE] * ((((*phaseField)(x1 + 2, x2 + 2, x3) - (*phaseField)(x1 - 2, x2 - 2, x3)) + ((*phaseField)(x1 + 2, x2 - 2, x3) - (*phaseField)(x1 - 2, x2 + 2, x3)))
+ // + (((*phaseField)(x1 + 2, x2, x3 - 2) - (*phaseField)(x1 - 2, x2, x3 + 2)) + ((*phaseField)(x1 + 2, x2, x3 + 2) - (*phaseField)(x1 - 2, x2, x3 - 2)))))
+ // + WEIGTH[N] * ((*phaseField)(x1 + 2, x2, x3) - (*phaseField)(x1 - 2, x2, x3))
+ // );
+ ////if (dX1_phi != NdX1_phi) {std::cout<<dX1_phi<<" "<< NdX1_phi<<std::endl;}
+
+ //LBMReal dX2_phi2 = 1.5 * ((
+ // WEIGTH[TNE] * ((((*phaseField)(x1 + 2, x2 + 2, x3 + 2) - (*phaseField)(x1 - 2, x2 - 2, x3 - 2)) + ((*phaseField)(x1 - 2, x2 + 2, x3 + 2) - (*phaseField)(x1 + 2, x2 - 2, x3 - 2)))
+ // + (((*phaseField)(x1 - 2, x2 + 2, x3 - 2) - (*phaseField)(x1 + 2, x2 - 2, x3 + 2)) + ((*phaseField)(x1 + 2, x2 + 2, x3 - 2) - (*phaseField)(x1 - 2, x2 - 2, x3 + 2))))
+ // + WEIGTH[NE] * ((((*phaseField)(x1 + 2, x2 + 2, x3) - (*phaseField)(x1 - 2, x2 - 2, x3)) + ((*phaseField)(x1 - 2, x2 + 2, x3) - (*phaseField)(x1 + 2, x2 - 2, x3)))
+ // + (((*phaseField)(x1, x2 + 2, x3 - 2) - (*phaseField)(x1, x2 - 2, x3 + 2)) + ((*phaseField)(x1, x2 + 2, x3 + 2) - (*phaseField)(x1, x2 - 2, x3 - 2)))))
+ // + WEIGTH[N] * ((*phaseField)(x1, x2 + 2, x3) - (*phaseField)(x1, x2 - 2, x3))
+ // );
+
+ //LBMReal dX3_phi2 = 1.5 * ((
+ // WEIGTH[TNE] * ((((*phaseField)(x1 + 2, x2 + 2, x3 + 2) - (*phaseField)(x1 - 2, x2 - 2, x3 - 2)) + ((*phaseField)(x1 - 2, x2 + 2, x3 + 2) - (*phaseField)(x1 + 2, x2 - 2, x3 - 2)))
+ // + (((*phaseField)(x1 - 2, x2 - 2, x3 + 2) - (*phaseField)(x1 + 2, x2 + 2, x3 - 2)) + ((*phaseField)(x1 + 2, x2 - 2, x3 + 2) - (*phaseField)(x1 - 2, x2 + 2, x3 - 2))))
+ // + WEIGTH[NE] * ((((*phaseField)(x1 + 2, x2, x3 + 2) - (*phaseField)(x1 - 2, x2, x3 - 2)) + ((*phaseField)(x1 - 2, x2, x3 + 2) - (*phaseField)(x1 + 2, x2, x3 - 2)))
+ // + (((*phaseField)(x1, x2 - 2, x3 + 2) - (*phaseField)(x1, x2 + 2, x3 - 2)) + ((*phaseField)(x1, x2 + 2, x3 + 2) - (*phaseField)(x1, x2 - 2, x3 - 2)))))
+ // + WEIGTH[N] * ((*phaseField)(x1, x2, x3 + 2) - (*phaseField)(x1, x2, x3 - 2))
+ // );
+
+ //dX1_phi = (2*dX1_phi -1*dX1_phi2);// 2 * dX1_phi - dX1_phi2;
+ //dX2_phi = (2*dX2_phi -1*dX2_phi2);// 2 * dX2_phi - dX2_phi2;
+ //dX3_phi = (2*dX3_phi -1*dX3_phi2);// 2 * dX3_phi - dX3_phi2;
+
LBMReal denom = sqrt(dX1_phi * dX1_phi + dX2_phi * dX2_phi + dX3_phi * dX3_phi) + 1e-9;
LBMReal normX1 = dX1_phi/denom;
@@ -312,20 +439,20 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
//----------- Calculating Macroscopic Values -------------
LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
- 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();
+ 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);
- LBMReal rho_m = 1.0 / densityRatio;
- forcingX1 = forcingX1 * (rho - rho_m);
- forcingX2 = forcingX2 * (rho - rho_m);
- forcingX3 = forcingX3 * (rho - rho_m);
+ forcingX1 = muForcingX1.Eval();
+ forcingX2 = muForcingX2.Eval();
+ forcingX3 = muForcingX3.Eval();
+ LBMReal rho_m = 1.0 / densityRatio;
+ forcingX1 = forcingX1 * (rho - rho_m);
+ forcingX2 = forcingX2 * (rho - rho_m);
+ forcingX3 = forcingX3 * (rho - rho_m);
+ }
////Incompressible Kernal
mfbbc = (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3)/rho*c3;
@@ -387,69 +514,272 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
vvz += mu * dX3_phi * c1o2;
+
+ ////Velocity filter 14.04.2021
+ // LBMReal lap_vx, lap_vy,lap_vz;
+ // {
+ // LBMReal sum = 0.0;
+ // sum += WEIGTH[TNE] * (((((*velocityX)(x1+1, x2+1, x3+1) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1 - 1, x2 - 1, x3 - 1) - (*velocityX)(x1, x2, x3))) + (((*velocityX)(x1 + 1, x2 + 1, x3 - 1) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1 + 1, x2 - 1, x3 + 1) - (*velocityX)(x1, x2, x3))))
+ // + ((((*velocityX)(x1 + 1, x2 - 1, x3 + 1) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1 - 1, x2 + 1, x3 - 1) - (*velocityX)(x1, x2, x3))) + (((*velocityX)(x1 - 1, x2 + 1, x3 + 1) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1 + 1, x2 - 1, x3 - 1) - (*velocityX)(x1, x2, x3)))));
+ // sum += WEIGTH[TN] * (
+ // ((((*velocityX)(x1 + 1, x2 + 1, x3 ) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1 - 1, x2 - 1, x3) - (*velocityX)(x1, x2, x3))) + (((*velocityX)(x1 + 1, x2 - 1, x3) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1 - 1, x2 + 1, x3) - (*velocityX)(x1, x2, x3))))
+ // + ((((*velocityX)(x1 + 1, x2 , x3+1) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1 - 1, x2 , x3-1) - (*velocityX)(x1, x2, x3))) + (((*velocityX)(x1 +1 , x2 , x3-1) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1 - 1, x2, x3 + 1) - (*velocityX)(x1, x2, x3))))
+ // + ((((*velocityX)(x1 , x2+1, x3 + 1) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1, x2 - 1, x3 - 1) - (*velocityX)(x1, x2, x3))) + (((*velocityX)(x1, x2 + 1, x3 - 1) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1, x2 - 1, x3 + 1) - (*velocityX)(x1, x2, x3))))
+ // );
+ // sum += WEIGTH[T] * (
+ // (((*velocityX)(x1-1, x2 , x3 ) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1 + 1, x2, x3) - (*velocityX)(x1, x2, x3)))
+ // + (((*velocityX)(x1 , x2-1, x3) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1, x2 + 1, x3) - (*velocityX)(x1, x2, x3)))
+ // + (((*velocityX)(x1, x2, x3-1) - (*velocityX)(x1, x2, x3)) + ((*velocityX)(x1, x2, x3+1) - (*velocityX)(x1, x2, x3)))
+ // );
+ // //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+ // // sum += WEIGTH[k] * (phi[k] - phi[REST]);
+ // //}
+ // lap_vx=6.0 * sum;
+
+ //sum = 0.0;
+ //sum += WEIGTH[TNE] * (((((*velocityY)(x1 + 1, x2 + 1, x3 + 1) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1 - 1, x2 - 1, x3 - 1) - (*velocityY)(x1, x2, x3))) + (((*velocityY)(x1 + 1, x2 + 1, x3 - 1) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1 + 1, x2 - 1, x3 + 1) - (*velocityY)(x1, x2, x3))))
+ // + ((((*velocityY)(x1 + 1, x2 - 1, x3 + 1) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1 - 1, x2 + 1, x3 - 1) - (*velocityY)(x1, x2, x3))) + (((*velocityY)(x1 - 1, x2 + 1, x3 + 1) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1 + 1, x2 - 1, x3 - 1) - (*velocityY)(x1, x2, x3)))));
+ //sum += WEIGTH[TN] * (
+ // ((((*velocityY)(x1 + 1, x2 + 1, x3) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1 - 1, x2 - 1, x3) - (*velocityY)(x1, x2, x3))) + (((*velocityY)(x1 + 1, x2 - 1, x3) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1 - 1, x2 + 1, x3) - (*velocityY)(x1, x2, x3))))
+ // + ((((*velocityY)(x1 + 1, x2, x3 + 1) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1 - 1, x2, x3 - 1) - (*velocityY)(x1, x2, x3))) + (((*velocityY)(x1 + 1, x2, x3 - 1) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1 - 1, x2, x3 + 1) - (*velocityY)(x1, x2, x3))))
+ // + ((((*velocityY)(x1, x2 + 1, x3 + 1) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1, x2 - 1, x3 - 1) - (*velocityY)(x1, x2, x3))) + (((*velocityY)(x1, x2 + 1, x3 - 1) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1, x2 - 1, x3 + 1) - (*velocityY)(x1, x2, x3))))
+ // );
+ //sum += WEIGTH[T] * (
+ // (((*velocityY)(x1 - 1, x2, x3) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1 + 1, x2, x3) - (*velocityY)(x1, x2, x3)))
+ // + (((*velocityY)(x1, x2 - 1, x3) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1, x2 + 1, x3) - (*velocityY)(x1, x2, x3)))
+ // + (((*velocityY)(x1, x2, x3 - 1) - (*velocityY)(x1, x2, x3)) + ((*velocityY)(x1, x2, x3 + 1) - (*velocityY)(x1, x2, x3)))
+ // );
+
+ //lap_vy = 6.0 * sum;
+
+ //sum = 0.0;
+ //sum += WEIGTH[TNE] * (((((*velocityZ)(x1 + 1, x2 + 1, x3 + 1) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1 - 1, x2 - 1, x3 - 1) - (*velocityZ)(x1, x2, x3))) + (((*velocityZ)(x1 + 1, x2 + 1, x3 - 1) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1 + 1, x2 - 1, x3 + 1) - (*velocityZ)(x1, x2, x3))))
+ // + ((((*velocityZ)(x1 + 1, x2 - 1, x3 + 1) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1 - 1, x2 + 1, x3 - 1) - (*velocityZ)(x1, x2, x3))) + (((*velocityZ)(x1 - 1, x2 + 1, x3 + 1) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1 + 1, x2 - 1, x3 - 1) - (*velocityZ)(x1, x2, x3)))));
+ //sum += WEIGTH[TN] * (
+ // ((((*velocityZ)(x1 + 1, x2 + 1, x3) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1 - 1, x2 - 1, x3) - (*velocityZ)(x1, x2, x3))) + (((*velocityZ)(x1 + 1, x2 - 1, x3) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1 - 1, x2 + 1, x3) - (*velocityZ)(x1, x2, x3))))
+ // + ((((*velocityZ)(x1 + 1, x2, x3 + 1) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1 - 1, x2, x3 - 1) - (*velocityZ)(x1, x2, x3))) + (((*velocityZ)(x1 + 1, x2, x3 - 1) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1 - 1, x2, x3 + 1) - (*velocityZ)(x1, x2, x3))))
+ // + ((((*velocityZ)(x1, x2 + 1, x3 + 1) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1, x2 - 1, x3 - 1) - (*velocityZ)(x1, x2, x3))) + (((*velocityZ)(x1, x2 + 1, x3 - 1) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1, x2 - 1, x3 + 1) - (*velocityZ)(x1, x2, x3))))
+ // );
+ //sum += WEIGTH[T] * (
+ // (((*velocityZ)(x1 - 1, x2, x3) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1 + 1, x2, x3) - (*velocityZ)(x1, x2, x3)))
+ // + (((*velocityZ)(x1, x2 - 1, x3) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1, x2 + 1, x3) - (*velocityZ)(x1, x2, x3)))
+ // + (((*velocityZ)(x1, x2, x3 - 1) - (*velocityZ)(x1, x2, x3)) + ((*velocityZ)(x1, x2, x3 + 1) - (*velocityZ)(x1, x2, x3)))
+ // );
+
+ //lap_vz = 6.0 * sum;
+
+ // }
+
+ // if (lap_vx != 0.0) {
+ // lap_vx = lap_vx;
+ // }
+
///----Classic source term 8.4.2021
+ LBMReal vvxF, vvyF, vvzF;
+ vvxF = vvx;//-2*c1o24 * lap_vx;//
+ vvyF = vvy;//-2*c1o24 * lap_vy;//
+ vvzF = vvz;//-2*c1o24 * lap_vz;//
+
+// vvxF = 1.2* vvx- 0.2*0.5 * ((*velocityX)(x1 - 1, x2, x3) + (*velocityX)(x1 + 1, x2, x3));
+// vvyF = 1.2 *vvy- 0.2*0.5* ((*velocityY)(x1 , x2-1, x3) + (*velocityY)(x1 , x2+1, x3));
+// vvzF = 1.2 *vvz-0.2*0.5* ((*velocityZ)(x1 , x2, x3-1) + (*velocityZ)(x1 , x2, x3+1));
+ //if (vvxF != vvx) {
+ // vvxF = vvxF;
+ //}
+ LBMReal weightGrad = 1.0-denom*denom/(denom*denom+0.0001*0.001);
+ LBMReal dX1_phiF = dX1_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX1;
+ LBMReal dX2_phiF = dX2_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX2;
+ LBMReal dX3_phiF = dX3_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX3;
+
+ //dX1_phiF *= 1.2;
+ //dX2_phiF *= 1.2;
+ //dX3_phiF *= 1.2;
+
+ //LBMReal gradFD = sqrt(dX1_phi * dX1_phi + dX2_phi * dX2_phi + dX3_phi * dX3_phi);
+ //LBMReal gradPhi = (1.0 - phi[REST]) * (phi[REST]);
+ //gradPhi = (gradPhi > gradFD) ? gradPhi : gradFD;
+ //dX1_phiF = gradPhi * normX1;
+ // dX2_phiF = gradPhi * normX2;
+ // dX3_phiF = gradPhi * normX3;
+
LBMReal ux2;
LBMReal uy2;
LBMReal uz2;
- ux2 = vvx * vvx;
- uy2 = vvy * vvy;
- uz2 = vvz * vvz;
+ ux2 = vvxF * vvxF;
+ uy2 = vvyF * vvyF;
+ uz2 = vvzF * vvzF;
LBMReal forcingTerm[D3Q27System::ENDF + 1];
for (int dir = STARTF; dir <= (FENDDIR); dir++) {
- LBMReal velProd = DX1[dir] * vvx + DX2[dir] * vvy + DX3[dir] * vvz;
+ LBMReal velProd = DX1[dir] * vvxF + DX2[dir] * vvyF + DX3[dir] * vvzF;
LBMReal velSq1 = velProd * velProd;
- LBMReal gamma = WEIGTH[dir] * (1.0 + 3 * velProd + 4.5 * velSq1 - 1.5 * (ux2 + uy2 + uz2));
+ LBMReal gamma = WEIGTH[dir] * (1.0 + 3 * velProd + (4.5 * velSq1 - 1.5 * (ux2 + uy2 + uz2)));
LBMReal fac1 = (gamma - WEIGTH[dir]) * c1o3 * rhoToPhi;
forcingTerm[dir] =
- (-vvx) * (fac1 * dX1_phi ) +
- (-vvy) * (fac1 * dX2_phi ) +
- (-vvz) * (fac1 * dX3_phi ) +
- (DX1[dir]) * (fac1 * dX1_phi ) +
- (DX2[dir]) * (fac1 * dX2_phi ) +
- (DX3[dir]) * (fac1 * dX3_phi );
+ (-vvxF) * (fac1 * dX1_phiF ) +
+ (-vvyF) * (fac1 * dX2_phiF ) +
+ (-vvzF) * (fac1 * dX3_phiF ) +
+ (DX1[dir]) * (fac1 * dX1_phiF ) +
+ (DX2[dir]) * (fac1 * dX2_phiF ) +
+ (DX3[dir]) * (fac1 * dX3_phiF );
+
+ //LBMReal biDif= (-((*phaseField)(x1 + 2 * DX1[dir], x2 + 2 * DX2[dir], x3 + 2 * DX3[dir])) + 4 * ((*phaseField)(x1 + DX1[dir], x2 + DX2[dir], x3 + DX3[dir]))
+ // - 3*((*phaseField)(x1 , x2 , x3 )) )*0.5;
+ //LBMReal ceDif = (((*phaseField)(x1 + DX1[dir], x2 + DX2[dir], x3 + DX3[dir])) - ((*phaseField)(x1 - DX1[dir], x2 - DX2[dir], x3 - DX3[dir]))) * 0.5;
+
+ ////ceDif = ((((*phaseField)(x1 + 2*DX1[dir], x2 + 2*DX2[dir], x3 + 2*DX3[dir])) - ((*phaseField)(x1 , x2 , x3 ))) * biDif < 0) ?
+ //// (!bcArray->isSolid(x1+2*DX1[dir], x2+2*DX2[dir], x3+2*DX3[dir]) && !bcArray->isUndefined(x1 + 2 * DX1[dir], x2 + 2 * DX2[dir], x3 + 2 * DX3[dir]) && !bcArray->isSolid(x1 + DX1[dir], x2 + DX2[dir], x3 + DX3[dir]) && !bcArray->isUndefined(x1 + DX1[dir], x2 + DX2[dir], x3 + DX3[dir]) && !bcArray->isSolid(x1 - DX1[dir], x2 - DX2[dir], x3 - DX3[dir]) && !bcArray->isUndefined(x1 - DX1[dir], x2 - DX2[dir], x3 - DX3[dir])) ?
+ //// (biDif+ceDif)*0.5 : ceDif: ceDif;
+
+ //ceDif = ((((*phaseField)(x1 + 2 * DX1[dir], x2 + 2 * DX2[dir], x3 + 2 * DX3[dir])) - ((*phaseField)(x1, x2, x3))) * biDif < 0) ? biDif : ceDif;
+
+ //forcingTerm[dir] =
+ // (-vvxF) * (fac1 * dX1_phiF) +
+ // (-vvyF) * (fac1 * dX2_phiF) +
+ // (-vvzF) * (fac1 * dX3_phiF) +
+ // fac1 * ceDif;//(((*phaseField)(x1 + DX1[dir], x2 + DX2[dir], x3 + DX3[dir])) - ((*phaseField)(x1 - DX1[dir], x2 - DX2[dir], x3 - DX3[dir]))) * 0.5;
+ // //( -((*phaseField)(x1 +2* DX1[dir], x2 + 2 * DX2[dir], x3 + 2 * DX3[dir])) + 5*((*phaseField)(x1 + DX1[dir], x2 + DX2[dir], x3 + DX3[dir]))
+ // //- 3*((*phaseField)(x1 , x2 , x3 )) - ((*phaseField)(x1 - DX1[dir], x2 - DX2[dir], x3 - DX3[dir])) )*0.25;
+
+
}
LBMReal gamma = WEIGTH[REST] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
LBMReal fac1 = (gamma - WEIGTH[REST]) * c1o3 * rhoToPhi;
- forcingTerm[REST] = (-vvx) * (fac1 * dX1_phi ) +
- (-vvy) * (fac1 * dX2_phi ) +
- (-vvz) * (fac1 * dX3_phi );
-
- mfcbb += 3.0 * ( 0.5 * forcingTerm[E]) / rho; //-(3.0*p1 - rho)*WEIGTH[E ];
- mfbcb += 3.0 * ( 0.5 * forcingTerm[N]) / rho; //-(3.0*p1 - rho)*WEIGTH[N ];
- mfbbc += 3.0 * ( 0.5 * forcingTerm[T]) / rho; //-(3.0*p1 - rho)*WEIGTH[T ];
- mfccb += 3.0 * ( 0.5 * forcingTerm[NE]) / rho; //-(3.0*p1 - rho)*WEIGTH[NE ];
- mfacb += 3.0 * ( 0.5 * forcingTerm[NW]) / rho; //-(3.0*p1 - rho)*WEIGTH[NW ];
- mfcbc += 3.0 * ( 0.5 * forcingTerm[TE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TE ];
- mfabc += 3.0 * ( 0.5 * forcingTerm[TW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TW ];
- mfbcc += 3.0 * ( 0.5 * forcingTerm[TN]) / rho; //-(3.0*p1 - rho)*WEIGTH[TN ];
- mfbac += 3.0 * ( 0.5 * forcingTerm[TS]) / rho; //-(3.0*p1 - rho)*WEIGTH[TS ];
- mfccc += 3.0 * ( 0.5 * forcingTerm[TNE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TNE];
- mfacc += 3.0 * ( 0.5 * forcingTerm[TNW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TNW];
- mfcac += 3.0 * ( 0.5 * forcingTerm[TSE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TSE];
- mfaac += 3.0 * ( 0.5 * forcingTerm[TSW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TSW];
- mfabb += 3.0 * ( 0.5 * forcingTerm[W]) / rho; //-(3.0*p1 - rho)*WEIGTH[W ];
- mfbab += 3.0 * ( 0.5 * forcingTerm[S]) / rho; //-(3.0*p1 - rho)*WEIGTH[S ];
- mfbba += 3.0 * ( 0.5 * forcingTerm[B]) / rho; //-(3.0*p1 - rho)*WEIGTH[B ];
- mfaab += 3.0 * ( 0.5 * forcingTerm[SW]) / rho; //-(3.0*p1 - rho)*WEIGTH[SW ];
- mfcab += 3.0 * ( 0.5 * forcingTerm[SE]) / rho; //-(3.0*p1 - rho)*WEIGTH[SE ];
- mfaba += 3.0 * ( 0.5 * forcingTerm[BW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BW ];
- mfcba += 3.0 * ( 0.5 * forcingTerm[BE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BE ];
- mfbaa += 3.0 * ( 0.5 * forcingTerm[BS]) / rho; //-(3.0*p1 - rho)*WEIGTH[BS ];
- mfbca += 3.0 * ( 0.5 * forcingTerm[BN]) / rho; //-(3.0*p1 - rho)*WEIGTH[BN ];
- mfaaa += 3.0 * ( 0.5 * forcingTerm[BSW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BSW];
- mfcaa += 3.0 * ( 0.5 * forcingTerm[BSE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BSE];
- mfaca += 3.0 * ( 0.5 * forcingTerm[BNW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BNW];
- mfcca += 3.0 * ( 0.5 * forcingTerm[BNE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BNE];
- mfbbb += 3.0 * ( 0.5 * forcingTerm[REST]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
+ forcingTerm[REST] = (-vvxF) * (fac1 * dX1_phiF ) +
+ (-vvyF) * (fac1 * dX2_phiF ) +
+ (-vvzF) * (fac1 * dX3_phiF );
- //--------------------------------------------------------
+ ////////
+ // LBMReal divAfterSource=
+ //( mfcbb + 3.0 * (0.5 * forcingTerm[E]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF) *(vvyF) +(vvzF) *(vvzF)-1)+
+ //( mfbcb + 3.0 * (0.5 * forcingTerm[N]) / rho ) *((vvxF) *(vvxF) +(vvyF-1)*(vvyF-1)+(vvzF) *(vvzF)-1)+
+ //( mfbbc + 3.0 * (0.5 * forcingTerm[T]) / rho ) *((vvxF) *(vvxF) +(vvyF) *(vvyF) +(vvzF-1)*(vvzF-1)-1)+
+ //( mfccb + 3.0 * (0.5 * forcingTerm[NE]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF-1)*(vvyF-1)+(vvzF) *(vvzF)-1)+
+ //( mfacb + 3.0 * (0.5 * forcingTerm[NW]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF-1)*(vvyF-1)+(vvzF) *(vvzF)-1)+
+ //( mfcbc + 3.0 * (0.5 * forcingTerm[TE]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF) *(vvyF) +(vvzF-1)*(vvzF-1)-1)+
+ //( mfabc + 3.0 * (0.5 * forcingTerm[TW]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF) *(vvyF) +(vvzF-1)*(vvzF-1)-1)+
+ //( mfbcc + 3.0 * (0.5 * forcingTerm[TN]) / rho ) *((vvxF) *(vvxF) +(vvyF-1)*(vvyF-1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfbac + 3.0 * (0.5 * forcingTerm[TS]) / rho ) *((vvxF) *(vvxF) +(vvyF+1)*(vvyF+1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfccc + 3.0 * (0.5 * forcingTerm[TNE]) / rho) *((vvxF-1)*(vvxF-1)+(vvyF-1)*(vvyF-1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfacc + 3.0 * (0.5 * forcingTerm[TNW]) / rho) *((vvxF+1)*(vvxF+1)+(vvyF-1)*(vvyF-1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfcac + 3.0 * (0.5 * forcingTerm[TSE]) / rho) *((vvxF-1)*(vvxF-1)+(vvyF+1)*(vvyF+1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfaac + 3.0 * (0.5 * forcingTerm[TSW]) / rho) *((vvxF+1)*(vvxF+1)+(vvyF+1)*(vvyF+1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfabb + 3.0 * (0.5 * forcingTerm[W]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF) *(vvyF) +(vvzF) *(vvzF)-1)+
+ //( mfbab + 3.0 * (0.5 * forcingTerm[S]) / rho ) *((vvxF) *(vvxF) +(vvyF+1)*(vvyF+1)+(vvzF) *(vvzF)-1)+
+ //( mfbba + 3.0 * (0.5 * forcingTerm[B]) / rho ) *((vvxF) *(vvxF) +(vvyF) *(vvyF) +(vvzF+1)*(vvzF+1)-1)+
+ //( mfaab + 3.0 * (0.5 * forcingTerm[SW]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF+1)*(vvyF+1)+(vvzF) *(vvzF)-1)+
+ //( mfcab + 3.0 * (0.5 * forcingTerm[SE]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF+1)*(vvyF+1)+(vvzF) *(vvzF)-1)+
+ //( mfaba + 3.0 * (0.5 * forcingTerm[BW]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF) *(vvyF) +(vvzF+1)*(vvzF+1)-1)+
+ //( mfcba + 3.0 * (0.5 * forcingTerm[BE]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF) *(vvyF) +(vvzF+1)*(vvzF+1)-1)+
+ //( mfbaa + 3.0 * (0.5 * forcingTerm[BS]) / rho ) *((vvxF) *(vvxF) +(vvyF+1)*(vvyF+1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfbca + 3.0 * (0.5 * forcingTerm[BN]) / rho ) *((vvxF) *(vvxF) +(vvyF-1)*(vvyF-1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfaaa + 3.0 * (0.5 * forcingTerm[BSW]) / rho) *((vvxF+1)*(vvxF+1)+(vvyF+1)*(vvyF+1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfcaa + 3.0 * (0.5 * forcingTerm[BSE]) / rho) *((vvxF-1)*(vvxF-1)+(vvyF+1)*(vvyF+1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfaca + 3.0 * (0.5 * forcingTerm[BNW]) / rho) *((vvxF+1)*(vvxF+1)+(vvyF-1)*(vvyF-1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfcca + 3.0 * (0.5 * forcingTerm[BNE]) / rho) *((vvxF-1)*(vvxF-1)+(vvyF-1)*(vvyF-1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfbbb + 3.0 * (0.5 * forcingTerm[REST]) / rho)*((vvxF)*(vvxF)+(vvyF)*(vvyF)+(vvzF)*(vvzF)-1);
+
+ // LBMReal divBeforeSource =
+ // (mfcbb) * ((vvxF - 1) * (vvxF - 1) + (vvyF) * (vvyF)+(vvzF) * (vvzF)-1) +
+ // (mfbcb) * ((vvxF) * (vvxF)+(vvyF - 1) * (vvyF - 1) + (vvzF) * (vvzF)-1) +
+ // (mfbbc) * ((vvxF) * (vvxF)+(vvyF) * (vvyF)+(vvzF - 1) * (vvzF - 1)-1) +
+ // (mfccb) * ((vvxF - 1) * (vvxF - 1) + (vvyF - 1) * (vvyF - 1) + (vvzF) * (vvzF)-1) +
+ // (mfacb) * ((vvxF + 1) * (vvxF + 1) + (vvyF - 1) * (vvyF - 1) + (vvzF) * (vvzF)-1) +
+ // (mfcbc) * ((vvxF - 1) * (vvxF - 1) + (vvyF) * (vvyF)+(vvzF - 1) * (vvzF - 1)-1) +
+ // (mfabc) * ((vvxF + 1) * (vvxF + 1) + (vvyF) * (vvyF)+(vvzF - 1) * (vvzF - 1)-1) +
+ // (mfbcc) * ((vvxF) * (vvxF)+(vvyF - 1) * (vvyF - 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfbac) * ((vvxF) * (vvxF)+(vvyF + 1) * (vvyF + 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfccc) * ((vvxF - 1) * (vvxF - 1) + (vvyF - 1) * (vvyF - 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfacc) * ((vvxF + 1) * (vvxF + 1) + (vvyF - 1) * (vvyF - 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfcac) * ((vvxF - 1) * (vvxF - 1) + (vvyF + 1) * (vvyF + 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfaac) * ((vvxF + 1) * (vvxF + 1) + (vvyF + 1) * (vvyF + 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfabb) * ((vvxF + 1) * (vvxF + 1) + (vvyF) * (vvyF)+(vvzF) * (vvzF)-1) +
+ // (mfbab) * ((vvxF) * (vvxF)+(vvyF + 1) * (vvyF + 1) + (vvzF) * (vvzF)-1) +
+ // (mfbba) * ((vvxF) * (vvxF)+(vvyF) * (vvyF)+(vvzF + 1) * (vvzF + 1)-1) +
+ // (mfaab) * ((vvxF + 1) * (vvxF + 1) + (vvyF + 1) * (vvyF + 1) + (vvzF) * (vvzF)-1) +
+ // (mfcab) * ((vvxF - 1) * (vvxF - 1) + (vvyF + 1) * (vvyF + 1) + (vvzF) * (vvzF)-1) +
+ // (mfaba) * ((vvxF + 1) * (vvxF + 1) + (vvyF) * (vvyF)+(vvzF + 1) * (vvzF + 1)-1) +
+ // (mfcba) * ((vvxF - 1) * (vvxF - 1) + (vvyF) * (vvyF)+(vvzF + 1) * (vvzF + 1)-1) +
+ // (mfbaa) * ((vvxF) * (vvxF)+(vvyF + 1) * (vvyF + 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfbca) * ((vvxF) * (vvxF)+(vvyF - 1) * (vvyF - 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfaaa) * ((vvxF + 1) * (vvxF + 1) + (vvyF + 1) * (vvyF + 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfcaa) * ((vvxF - 1) * (vvxF - 1) + (vvyF + 1) * (vvyF + 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfaca) * ((vvxF + 1) * (vvxF + 1) + (vvyF - 1) * (vvyF - 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfcca) * ((vvxF - 1) * (vvxF - 1) + (vvyF - 1) * (vvyF - 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfbbb) * ((vvxF) * (vvxF)+(vvyF) * (vvyF)+(vvzF) * (vvzF)-1);
+ //if (divAfterSource - divBeforeSource != 0 && phi[REST]>0.0001 && phi[REST]<0.999) {
+ // std::cout << phi[REST]<<" "<< divAfterSource << " " << divBeforeSource <<" "<< divAfterSource/ divBeforeSource << std::endl;
+ //}
+
+ //if (fabs(divAfterSource - divBeforeSource)/(fabs(divAfterSource) + fabs(divBeforeSource)+1e-10) > 1e-5) {
+ // LBMReal scaleDiv =0.95+(1-0.95)* (divBeforeSource) / (divBeforeSource - divAfterSource);
+
+ // forcingTerm[E] *=scaleDiv;
+ // forcingTerm[N] *=scaleDiv;
+ // forcingTerm[T] *=scaleDiv;
+ // forcingTerm[NE] *=scaleDiv;
+ // forcingTerm[NW] *=scaleDiv;
+ // forcingTerm[TE] *=scaleDiv;
+ // forcingTerm[TW] *=scaleDiv;
+ // forcingTerm[TN] *=scaleDiv;
+ // forcingTerm[TS] *=scaleDiv;
+ // forcingTerm[TNE] *=scaleDiv;
+ // forcingTerm[TNW] *=scaleDiv;
+ // forcingTerm[TSE] *=scaleDiv;
+ // forcingTerm[TSW] *=scaleDiv;
+ // forcingTerm[W] *=scaleDiv;
+ // forcingTerm[S] *=scaleDiv;
+ // forcingTerm[B] *=scaleDiv;
+ // forcingTerm[SW] *=scaleDiv;
+ // forcingTerm[SE] *=scaleDiv;
+ // forcingTerm[BW] *=scaleDiv;
+ // forcingTerm[BE] *=scaleDiv;
+ // forcingTerm[BS] *=scaleDiv;
+ // forcingTerm[BN] *=scaleDiv;
+ // forcingTerm[BSW] *=scaleDiv;
+ // forcingTerm[BSE] *=scaleDiv;
+ // forcingTerm[BNW] *=scaleDiv;
+ // forcingTerm[BNE] *=scaleDiv;
+ // forcingTerm[REST] *=scaleDiv;
+ //}
+ ////////
+
+
+ mfcbb +=3.0 * ( 0.5 * forcingTerm[E]) / rho; //-(3.0*p1 - rho)*WEIGTH[E ];
+ mfbcb +=3.0 * ( 0.5 * forcingTerm[N]) / rho; //-(3.0*p1 - rho)*WEIGTH[N ];
+ mfbbc +=3.0 * ( 0.5 * forcingTerm[T]) / rho; //-(3.0*p1 - rho)*WEIGTH[T ];
+ mfccb +=3.0 * ( 0.5 * forcingTerm[NE]) / rho; //-(3.0*p1 - rho)*WEIGTH[NE ];
+ mfacb +=3.0 * ( 0.5 * forcingTerm[NW]) / rho; //-(3.0*p1 - rho)*WEIGTH[NW ];
+ mfcbc +=3.0 * ( 0.5 * forcingTerm[TE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TE ];
+ mfabc +=3.0 * ( 0.5 * forcingTerm[TW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TW ];
+ mfbcc +=3.0 * ( 0.5 * forcingTerm[TN]) / rho; //-(3.0*p1 - rho)*WEIGTH[TN ];
+ mfbac +=3.0 * ( 0.5 * forcingTerm[TS]) / rho; //-(3.0*p1 - rho)*WEIGTH[TS ];
+ mfccc +=3.0 * ( 0.5 * forcingTerm[TNE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TNE];
+ mfacc +=3.0 * ( 0.5 * forcingTerm[TNW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TNW];
+ mfcac +=3.0 * ( 0.5 * forcingTerm[TSE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TSE];
+ mfaac +=3.0 * ( 0.5 * forcingTerm[TSW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TSW];
+ mfabb +=3.0 * ( 0.5 * forcingTerm[W]) / rho; //-(3.0*p1 - rho)*WEIGTH[W ];
+ mfbab +=3.0 * ( 0.5 * forcingTerm[S]) / rho; //-(3.0*p1 - rho)*WEIGTH[S ];
+ mfbba +=3.0 * ( 0.5 * forcingTerm[B]) / rho; //-(3.0*p1 - rho)*WEIGTH[B ];
+ mfaab +=3.0 * ( 0.5 * forcingTerm[SW]) / rho; //-(3.0*p1 - rho)*WEIGTH[SW ];
+ mfcab +=3.0 * ( 0.5 * forcingTerm[SE]) / rho; //-(3.0*p1 - rho)*WEIGTH[SE ];
+ mfaba +=3.0 * ( 0.5 * forcingTerm[BW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BW ];
+ mfcba +=3.0 * ( 0.5 * forcingTerm[BE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BE ];
+ mfbaa +=3.0 * ( 0.5 * forcingTerm[BS]) / rho; //-(3.0*p1 - rho)*WEIGTH[BS ];
+ mfbca +=3.0 * ( 0.5 * forcingTerm[BN]) / rho; //-(3.0*p1 - rho)*WEIGTH[BN ];
+ mfaaa +=3.0 * ( 0.5 * forcingTerm[BSW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BSW];
+ mfcaa +=3.0 * ( 0.5 * forcingTerm[BSE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BSE];
+ mfaca +=3.0 * ( 0.5 * forcingTerm[BNW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BNW];
+ mfcca +=3.0 * ( 0.5 * forcingTerm[BNE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BNE];
+ mfbbb +=3.0 * ( 0.5 * forcingTerm[REST]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
+ //--------------------------------------------------------
+ //////////End classic source term
//forcing
///////////////////////////////////////////////////////////////////////////////////////////
if (withForcing)
@@ -753,11 +1083,11 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
/////fourth order parameters; here only for test. Move out of loop!
LBMReal OxyyPxzz = 8.0 * (collFactorM - 2.0) * (OxxPyyPzz * (3.0 * collFactorM - 1.0) - 5.0 * collFactorM) / (8.0 * (5.0 - 2.0 * collFactorM) * collFactorM + OxxPyyPzz * (8.0 + collFactorM * (9.0 * collFactorM - 26.0)));
- LBMReal OxyyMxzz = 8.0 * (collFactorM - 2.0) * (collFactorM + OxxPyyPzz * (3.0 * collFactorM - 7.0)) / (OxxPyyPzz * (56.0 - 42.0 * collFactorM + 9.0 * collFactorM * collFactorM) - 8.0 * collFactorM);
- LBMReal Oxyz = 24.0 * (collFactorM - 2.0) * (4.0 * collFactorM * collFactorM + collFactorM * OxxPyyPzz * (18.0 - 13.0 * collFactorM) + OxxPyyPzz * OxxPyyPzz * (2.0 + collFactorM * (6.0 * collFactorM - 11.0))) / (16.0 * collFactorM * collFactorM * (collFactorM - 6.0) - 2.0 * collFactorM * OxxPyyPzz * (216.0 + 5.0 * collFactorM * (9.0 * collFactorM - 46.0)) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (3.0 * collFactorM - 10.0) * (15.0 * collFactorM - 28.0) - 48.0));
- LBMReal A = (4.0 * collFactorM * collFactorM + 2.0 * collFactorM * OxxPyyPzz * (collFactorM - 6.0) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (10.0 - 3.0 * collFactorM) - 4.0)) / ((collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
+ LBMReal OxyyMxzz = 8.0 * (collFactorM - 2.0) * (collFactorM + OxxPyyPzz * (3.0 * collFactorM - 7.0)) / (OxxPyyPzz * (56.0 - 42.0 * collFactorM + 9.0 * collFactorM * collFactorM) - 8.0 * collFactorM);
+ LBMReal Oxyz = 24.0 * (collFactorM - 2.0) * (4.0 * collFactorM * collFactorM + collFactorM * OxxPyyPzz * (18.0 - 13.0 * collFactorM) + OxxPyyPzz * OxxPyyPzz * (2.0 + collFactorM * (6.0 * collFactorM - 11.0))) / (16.0 * collFactorM * collFactorM * (collFactorM - 6.0) - 2.0 * collFactorM * OxxPyyPzz * (216.0 + 5.0 * collFactorM * (9.0 * collFactorM - 46.0)) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (3.0 * collFactorM - 10.0) * (15.0 * collFactorM - 28.0) - 48.0));
+ LBMReal A = (4.0 * collFactorM * collFactorM + 2.0 * collFactorM * OxxPyyPzz * (collFactorM - 6.0) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (10.0 - 3.0 * collFactorM) - 4.0)) / ((collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
//FIXME: warning C4459: declaration of 'B' hides global declaration (message : see declaration of 'D3Q27System::B' )
- LBMReal BB = (4.0 * collFactorM * OxxPyyPzz * (9.0 * collFactorM - 16.0) - 4.0 * collFactorM * collFactorM - 2.0 * OxxPyyPzz * OxxPyyPzz * (2.0 + 9.0 * collFactorM * (collFactorM - 2.0))) / (3.0 * (collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
+ LBMReal BB = (4.0 * collFactorM * OxxPyyPzz * (9.0 * collFactorM - 16.0) - 4.0 * collFactorM * collFactorM - 2.0 * OxxPyyPzz * OxxPyyPzz * (2.0 + 9.0 * collFactorM * (collFactorM - 2.0))) / (3.0 * (collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
//Cum 4.
@@ -802,13 +1132,25 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
// mxxPyyPzz += c2o3 * rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz);
//17.03.2021 attempt for statililization by assymptotically vanishing bias
- LBMReal correctionScaling =0.0* rhoToPhi /rho;// +0.5;// (vx2 + vy2 + vz2) * 100;// +0.5;//(vx2 + vy2 + vz2)*1000;
- mxxPyyPzz += (1.0/3.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz)* correctionScaling; // As in Hesam's code
- mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy)* correctionScaling;
- mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
- mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
- mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
- mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
+ //LBMReal correctionScaling = rhoToPhi /rho;// +0.5;// (vx2 + vy2 + vz2) * 100;// +0.5;//(vx2 + vy2 + vz2)*1000;
+ //mxxPyyPzz += (1.0/3.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz)* correctionScaling; // As in Hesam's code
+ //mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy)* correctionScaling;
+ //mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
+ //mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
+ //mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
+ //mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
+
+
+ //14.04.2021 filtered velocity
+
+ //LBMReal correctionScaling = rhoToPhi / rho;// +0.5;// (vx2 + vy2 + vz2) * 100;// +0.5;//(vx2 + vy2 + vz2)*1000;
+ //mxxPyyPzz += (1.0 / 3.0) * (dX1_phi * vvxF + dX2_phi * vvyF + dX3_phi * vvzF) * correctionScaling; // As in Hesam's code
+ //mxxMyy += c1o3 * (dX1_phi * vvxF - dX2_phi * vvyF) * correctionScaling;
+ //mxxMzz += c1o3 * (dX1_phi * vvxF - dX3_phi * vvzF) * correctionScaling;
+ //mfabb += c1o6 * (dX2_phi * vvzF + dX3_phi * vvyF) * correctionScaling;
+ //mfbab += c1o6 * (dX1_phi * vvzF + dX3_phi * vvxF) * correctionScaling;
+ //mfbba += c1o6 * (dX1_phi * vvyF + dX2_phi * vvxF) * correctionScaling;
+
LBMReal dxux = -c1o2 * collFactorM * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz);
LBMReal dyuy = dxux + collFactorM * c3o2 * mxxMyy;
@@ -819,6 +1161,20 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
LBMReal Dyz = -three * collFactorM * mfabb;
////relax unfiltered
+ //! divergenceFilter 10.05.2021
+ LBMReal divMag= (1.0 - phi[REST]) * (phi[REST])*10*5*sqrt(fabs((OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz))));
+ // LBMReal divMag = 500 *500* 50*(fabs((OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz))))* (fabs((OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz))));
+ //LBMReal divMag = (dX1_phi * dxux) > 0 ? (dX1_phi * dxux) : 0;
+ //divMag += (dX2_phi * dyuy) > 0 ? (dX2_phi * dyuy) : 0;
+ //divMag += (dX3_phi * dzuz) > 0 ? (dX3_phi * dzuz) : 0;
+ //divMag *= 5000;
+ //divMag+= denom * 10 * 5 * sqrt(fabs((OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz))));
+ //LBMReal divMag = 5000 * (fabs(dX1_phi * dxux)+fabs(dX2_phi * dyuy)+fabs(dX3_phi * dzuz));
+ collFactorM = collFactorM / (1.0 + 3.0 * divMag);
+
+ collFactorM = (collFactorM > 1.0) ? collFactorM : 1.0;
+
+
mxxPyyPzz += OxxPyyPzz * (/*mfaaa*/ - mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
mxxMyy += collFactorM * (-mxxMyy) - 3. * (1. - c1o2 * collFactorM) * (vx2 * dxux - vy2 * dyuy);
mxxMzz += collFactorM * (-mxxMzz) - 3. * (1. - c1o2 * collFactorM) * (vx2 * dxux - vz2 * dzuz);
@@ -848,15 +1204,29 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
//applying phase field gradients second part:
//mxxPyyPzz += c2o3 * rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz);
- mxxPyyPzz += (1.0 / 3.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling; // As in Hesam's code
- mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy) * correctionScaling;
- mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
- mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
- mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
- mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
+ //mxxPyyPzz += (1.0 / 3.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling; // As in Hesam's code
+ //mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy) * correctionScaling;
+ //mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
+ //mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
+ //mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
+ //mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
+
+
+ //////updated pressure
+ //mfaaa += (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling;
+
+
+ //mxxPyyPzz += (1.0 / 3.0) * (dX1_phi * vvxF + dX2_phi * vvyF + dX3_phi * vvzF) * correctionScaling; // As in Hesam's code
+ //mxxMyy += c1o3 * (dX1_phi * vvxF - dX2_phi * vvyF) * correctionScaling;
+ //mxxMzz += c1o3 * (dX1_phi * vvxF - dX3_phi * vvzF) * correctionScaling;
+ //mfabb += c1o6 * (dX2_phi * vvzF + dX3_phi * vvyF) * correctionScaling;
+ //mfbab += c1o6 * (dX1_phi * vvzF + dX3_phi * vvxF) * correctionScaling;
+ //mfbba += c1o6 * (dX1_phi * vvyF + dX2_phi * vvxF) * correctionScaling;
+
+
+ //////updated pressure
+ //mfaaa += (dX1_phi * vvxF + dX2_phi * vvyF + dX3_phi * vvzF) * correctionScaling;
- ////updated pressure
- mfaaa += (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling;
mxxPyyPzz += mfaaa;//12.03.21 shifted by mfaaa
// mxxPyyPzz = mfaaa; //12.03.21 reguarized pressure !?
@@ -2542,7 +2912,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
}
}
}
- }
+
dataSet->setPhaseField(divU);
}
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp
index 102af6035..cc90f0fc9 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp
@@ -289,6 +289,7 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
//-1 0 1
findNeighbors(phaseField, x1, x2, x3);
+ findNeighbors2(phaseField2, x1, x2, x3);
LBMReal mfcbb = (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3);
LBMReal mfbcb = (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3);
@@ -328,10 +329,22 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
LBMReal dX2_phi = gradX2_phi();
LBMReal dX3_phi = gradX3_phi();
- LBMReal denom = sqrt(dX1_phi * dX1_phi + dX2_phi * dX2_phi + dX3_phi * dX3_phi) + 1e-9;
- LBMReal normX1 = dX1_phi/denom;
- LBMReal normX2 = dX2_phi/denom;
- LBMReal normX3 = dX3_phi/denom;
+ LBMReal dX1_phi2 = gradX1_phi2();
+ LBMReal dX2_phi2 = gradX2_phi2();
+ LBMReal dX3_phi2 = gradX3_phi2();
+
+
+ LBMReal denom2 = sqrt(dX1_phi * dX1_phi+ dX1_phi2 * dX1_phi2 + dX2_phi * dX2_phi + dX2_phi2 * dX2_phi2 + dX3_phi * dX3_phi+ dX3_phi2 * dX3_phi2) + 1e-9;
+ LBMReal normX1 = (dX1_phi-dX1_phi2)/denom2;
+ LBMReal normX2 = (dX2_phi-dX2_phi2)/denom2;
+ LBMReal normX3 = (dX3_phi-dX3_phi2)/denom2;
+
+ //LBMReal denom = sqrt(dX1_phi * dX1_phi + dX2_phi * dX2_phi + dX3_phi * dX3_phi) + 1e-9;
+ //LBMReal normX1 = dX1_phi / denom;
+ //LBMReal normX2 = dX2_phi / denom;
+ //LBMReal normX3 = dX3_phi / denom;
+
+
collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
@@ -341,51 +354,38 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
//----------- Calculating Macroscopic Values -------------
LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
- 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();
-
- LBMReal rho_m = 1.0 / densityRatio;
- forcingX1 = forcingX1 * (rho - rho_m);
- forcingX2 = forcingX2 * (rho - rho_m);
- forcingX3 = forcingX3 * (rho - rho_m);
-
////Incompressible Kernal
- mfbbc = (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3)/rho;
- mfbcb = (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) / rho;
- mfccb = (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3) / rho;
- mfacb = (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3) / rho;
- mfcbb = (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) / rho;
- mfcbc = (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3) / rho;
- mfabc = (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3) / rho;
- mfbcc = (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3) / rho;
- mfbac = (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3) / rho;
- mfccc = (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3) / rho;
- mfacc = (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3) / rho;
- mfcac = (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3) / rho;
- mfaac = (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3) / rho;
-
- mfabb = (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3) / rho;
- mfbab = (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3) / rho;
- mfbba = (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p) / rho;
- mfaab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3) / rho;
- mfcab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3) / rho;
- mfaba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p) / rho;
- mfcba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p) / rho;
- mfbaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p) / rho;
- mfbca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p) / rho;
- mfaaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p) / rho;
- mfcaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p) / rho;
- mfaca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) / rho;
- mfcca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) / rho;
-
- mfbbb = (*this->zeroDistributionsF)(x1, x2, x3) / rho;
+ mfbbc = (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3) / rho * c3;
+ mfbcb = (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) / rho * c3;
+ mfccb = (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3) / rho * c3;
+ mfacb = (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3) / rho * c3;
+ mfcbb = (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) / rho * c3;
+ mfcbc = (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3) / rho * c3;
+ mfabc = (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3) / rho * c3;
+ mfbcc = (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3) / rho * c3;
+ mfbac = (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3) / rho * c3;
+ mfccc = (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3) / rho * c3;
+ mfacc = (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3) / rho * c3;
+ mfcac = (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3) / rho * c3;
+ mfaac = (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3) / rho * c3;
+
+ mfabb = (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3) / rho * c3;
+ mfbab = (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3) / rho * c3;
+ mfbba = (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p) / rho * c3;
+ mfaab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3) / rho * c3;
+ mfcab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3) / rho * c3;
+ mfaba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p) / rho * c3;
+ mfcba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p) / rho * c3;
+ mfbaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p) / rho * c3;
+ mfbca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p) / rho * c3;
+ mfaaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p) / rho * c3;
+ mfcaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p) / rho * c3;
+ mfaca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) / rho * c3;
+ mfcca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) / rho * c3;
+
+ mfbbb = (*this->zeroDistributionsF)(x1, x2, x3) / rho * c3;
+
LBMReal m0, m1, m2;
LBMReal rhoRef=c1;
@@ -405,11 +405,240 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
(((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
(mfbbc - mfbba))/rhoRef;
+
+ 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();
+
+ //LBMReal rho_m = 1.0 / densityRatio;
+ //forcingX1 = forcingX1 * (rho - rho_m);
+ //forcingX2 = forcingX2 * (rho - rho_m);
+ //forcingX3 = forcingX3 * (rho - rho_m);
+ vvx += forcingX1 * deltaT * 0.5; // X
+ vvy += forcingX2 * deltaT * 0.5; // Y
+ vvz += forcingX3 * deltaT * 0.5; // Z
+
+ }
+
+
///surface tension force
vvx += mu * dX1_phi*c1o2;
- vvy += mu * dX2_phi * c1o2;
+ vvy += mu * dX2_phi * c1o2 ;
vvz += mu * dX3_phi * c1o2;
+ //////classic source term
+ ///----Classic source term 8.4.2021
+
+ LBMReal vvxF, vvyF, vvzF;
+ vvxF = vvx;//-2*c1o24 * lap_vx;//
+ vvyF = vvy;//-2*c1o24 * lap_vy;//
+ vvzF = vvz;//-2*c1o24 * lap_vz;//
+
+// vvxF = 1.2* vvx- 0.2*0.5 * ((*velocityX)(x1 - 1, x2, x3) + (*velocityX)(x1 + 1, x2, x3));
+// vvyF = 1.2 *vvy- 0.2*0.5* ((*velocityY)(x1 , x2-1, x3) + (*velocityY)(x1 , x2+1, x3));
+// vvzF = 1.2 *vvz-0.2*0.5* ((*velocityZ)(x1 , x2, x3-1) + (*velocityZ)(x1 , x2, x3+1));
+ //if (vvxF != vvx) {
+ // vvxF = vvxF;
+ //}
+ LBMReal weightGrad = 1.0;// -denom * denom / (denom * denom + 0.0001 * 0.001);
+ LBMReal dX1_phiF = dX1_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX1;
+ LBMReal dX2_phiF = dX2_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX2;
+ LBMReal dX3_phiF = dX3_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX3;
+
+ //dX1_phiF *= 1.2;
+ //dX2_phiF *= 1.2;
+ //dX3_phiF *= 1.2;
+
+ //LBMReal gradFD = sqrt(dX1_phi * dX1_phi + dX2_phi * dX2_phi + dX3_phi * dX3_phi);
+ //LBMReal gradPhi = (1.0 - phi[REST]) * (phi[REST]);
+ //gradPhi = (gradPhi > gradFD) ? gradPhi : gradFD;
+ //dX1_phiF = gradPhi * normX1;
+ // dX2_phiF = gradPhi * normX2;
+ // dX3_phiF = gradPhi * normX3;
+
+ LBMReal ux2;
+ LBMReal uy2;
+ LBMReal uz2;
+ ux2 = vvxF * vvxF;
+ uy2 = vvyF * vvyF;
+ uz2 = vvzF * vvzF;
+ LBMReal forcingTerm[D3Q27System::ENDF + 1];
+ for (int dir = STARTF; dir <= (FENDDIR); dir++) {
+ LBMReal velProd = DX1[dir] * vvxF + DX2[dir] * vvyF + DX3[dir] * vvzF;
+ LBMReal velSq1 = velProd * velProd;
+ LBMReal gamma = WEIGTH[dir] * (1.0 + 3 * velProd + (4.5 * velSq1 - 1.5 * (ux2 + uy2 + uz2)));
+
+ //LBMReal fac1 = (gamma - WEIGTH[dir]) * c1o3 * rhoToPhi;
+
+ //forcingTerm[dir] =
+ // (-vvxF) * (fac1 * dX1_phiF) +
+ // (-vvyF) * (fac1 * dX2_phiF) +
+ // (-vvzF) * (fac1 * dX3_phiF) +
+ // (DX1[dir]) * (fac1 * dX1_phiF) +
+ // (DX2[dir]) * (fac1 * dX2_phiF) +
+ // (DX3[dir]) * (fac1 * dX3_phiF);
+
+
+ LBMReal fac1 = (gamma - WEIGTH[dir]) * c1o3 ;
+
+ forcingTerm[dir] =
+ (-vvxF) * (fac1 * (dX1_phiF * rhoH + dX2_phi2 * rhoL)) +
+ (-vvyF) * (fac1 * (dX2_phiF * rhoH + dX2_phi2 * rhoL)) +
+ (-vvzF) * (fac1 * (dX3_phiF * rhoH + dX3_phi2 * rhoL)) +
+ (DX1[dir]) * (fac1 * (dX1_phiF * rhoH + dX2_phi2 * rhoL)) +
+ (DX2[dir]) * (fac1 * (dX2_phiF * rhoH + dX2_phi2 * rhoL)) +
+ (DX3[dir]) * (fac1 * (dX3_phiF * rhoH + dX3_phi2 * rhoL));
+
+
+
+ }
+
+ LBMReal gamma = WEIGTH[REST] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
+ LBMReal fac1 = (gamma - WEIGTH[REST]) * c1o3 * rhoToPhi;
+ forcingTerm[REST] = (-vvxF) * (fac1 * (dX1_phiF * rhoH + dX2_phi2 * rhoL)) +
+ (-vvyF) * (fac1 * (dX2_phiF * rhoH + dX2_phi2 * rhoL)) +
+ (-vvzF) * (fac1 * (dX3_phiF * rhoH + dX3_phi2 * rhoL));
+
+ ////////
+ // LBMReal divAfterSource=
+ //( mfcbb + 3.0 * (0.5 * forcingTerm[E]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF) *(vvyF) +(vvzF) *(vvzF)-1)+
+ //( mfbcb + 3.0 * (0.5 * forcingTerm[N]) / rho ) *((vvxF) *(vvxF) +(vvyF-1)*(vvyF-1)+(vvzF) *(vvzF)-1)+
+ //( mfbbc + 3.0 * (0.5 * forcingTerm[T]) / rho ) *((vvxF) *(vvxF) +(vvyF) *(vvyF) +(vvzF-1)*(vvzF-1)-1)+
+ //( mfccb + 3.0 * (0.5 * forcingTerm[NE]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF-1)*(vvyF-1)+(vvzF) *(vvzF)-1)+
+ //( mfacb + 3.0 * (0.5 * forcingTerm[NW]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF-1)*(vvyF-1)+(vvzF) *(vvzF)-1)+
+ //( mfcbc + 3.0 * (0.5 * forcingTerm[TE]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF) *(vvyF) +(vvzF-1)*(vvzF-1)-1)+
+ //( mfabc + 3.0 * (0.5 * forcingTerm[TW]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF) *(vvyF) +(vvzF-1)*(vvzF-1)-1)+
+ //( mfbcc + 3.0 * (0.5 * forcingTerm[TN]) / rho ) *((vvxF) *(vvxF) +(vvyF-1)*(vvyF-1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfbac + 3.0 * (0.5 * forcingTerm[TS]) / rho ) *((vvxF) *(vvxF) +(vvyF+1)*(vvyF+1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfccc + 3.0 * (0.5 * forcingTerm[TNE]) / rho) *((vvxF-1)*(vvxF-1)+(vvyF-1)*(vvyF-1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfacc + 3.0 * (0.5 * forcingTerm[TNW]) / rho) *((vvxF+1)*(vvxF+1)+(vvyF-1)*(vvyF-1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfcac + 3.0 * (0.5 * forcingTerm[TSE]) / rho) *((vvxF-1)*(vvxF-1)+(vvyF+1)*(vvyF+1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfaac + 3.0 * (0.5 * forcingTerm[TSW]) / rho) *((vvxF+1)*(vvxF+1)+(vvyF+1)*(vvyF+1)+(vvzF-1)*(vvzF-1)-1)+
+ //( mfabb + 3.0 * (0.5 * forcingTerm[W]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF) *(vvyF) +(vvzF) *(vvzF)-1)+
+ //( mfbab + 3.0 * (0.5 * forcingTerm[S]) / rho ) *((vvxF) *(vvxF) +(vvyF+1)*(vvyF+1)+(vvzF) *(vvzF)-1)+
+ //( mfbba + 3.0 * (0.5 * forcingTerm[B]) / rho ) *((vvxF) *(vvxF) +(vvyF) *(vvyF) +(vvzF+1)*(vvzF+1)-1)+
+ //( mfaab + 3.0 * (0.5 * forcingTerm[SW]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF+1)*(vvyF+1)+(vvzF) *(vvzF)-1)+
+ //( mfcab + 3.0 * (0.5 * forcingTerm[SE]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF+1)*(vvyF+1)+(vvzF) *(vvzF)-1)+
+ //( mfaba + 3.0 * (0.5 * forcingTerm[BW]) / rho ) *((vvxF+1)*(vvxF+1)+(vvyF) *(vvyF) +(vvzF+1)*(vvzF+1)-1)+
+ //( mfcba + 3.0 * (0.5 * forcingTerm[BE]) / rho ) *((vvxF-1)*(vvxF-1)+(vvyF) *(vvyF) +(vvzF+1)*(vvzF+1)-1)+
+ //( mfbaa + 3.0 * (0.5 * forcingTerm[BS]) / rho ) *((vvxF) *(vvxF) +(vvyF+1)*(vvyF+1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfbca + 3.0 * (0.5 * forcingTerm[BN]) / rho ) *((vvxF) *(vvxF) +(vvyF-1)*(vvyF-1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfaaa + 3.0 * (0.5 * forcingTerm[BSW]) / rho) *((vvxF+1)*(vvxF+1)+(vvyF+1)*(vvyF+1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfcaa + 3.0 * (0.5 * forcingTerm[BSE]) / rho) *((vvxF-1)*(vvxF-1)+(vvyF+1)*(vvyF+1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfaca + 3.0 * (0.5 * forcingTerm[BNW]) / rho) *((vvxF+1)*(vvxF+1)+(vvyF-1)*(vvyF-1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfcca + 3.0 * (0.5 * forcingTerm[BNE]) / rho) *((vvxF-1)*(vvxF-1)+(vvyF-1)*(vvyF-1)+(vvzF+1)*(vvzF+1)-1)+
+ //( mfbbb + 3.0 * (0.5 * forcingTerm[REST]) / rho)*((vvxF)*(vvxF)+(vvyF)*(vvyF)+(vvzF)*(vvzF)-1);
+
+ // LBMReal divBeforeSource =
+ // (mfcbb) * ((vvxF - 1) * (vvxF - 1) + (vvyF) * (vvyF)+(vvzF) * (vvzF)-1) +
+ // (mfbcb) * ((vvxF) * (vvxF)+(vvyF - 1) * (vvyF - 1) + (vvzF) * (vvzF)-1) +
+ // (mfbbc) * ((vvxF) * (vvxF)+(vvyF) * (vvyF)+(vvzF - 1) * (vvzF - 1)-1) +
+ // (mfccb) * ((vvxF - 1) * (vvxF - 1) + (vvyF - 1) * (vvyF - 1) + (vvzF) * (vvzF)-1) +
+ // (mfacb) * ((vvxF + 1) * (vvxF + 1) + (vvyF - 1) * (vvyF - 1) + (vvzF) * (vvzF)-1) +
+ // (mfcbc) * ((vvxF - 1) * (vvxF - 1) + (vvyF) * (vvyF)+(vvzF - 1) * (vvzF - 1)-1) +
+ // (mfabc) * ((vvxF + 1) * (vvxF + 1) + (vvyF) * (vvyF)+(vvzF - 1) * (vvzF - 1)-1) +
+ // (mfbcc) * ((vvxF) * (vvxF)+(vvyF - 1) * (vvyF - 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfbac) * ((vvxF) * (vvxF)+(vvyF + 1) * (vvyF + 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfccc) * ((vvxF - 1) * (vvxF - 1) + (vvyF - 1) * (vvyF - 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfacc) * ((vvxF + 1) * (vvxF + 1) + (vvyF - 1) * (vvyF - 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfcac) * ((vvxF - 1) * (vvxF - 1) + (vvyF + 1) * (vvyF + 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfaac) * ((vvxF + 1) * (vvxF + 1) + (vvyF + 1) * (vvyF + 1) + (vvzF - 1) * (vvzF - 1)-1) +
+ // (mfabb) * ((vvxF + 1) * (vvxF + 1) + (vvyF) * (vvyF)+(vvzF) * (vvzF)-1) +
+ // (mfbab) * ((vvxF) * (vvxF)+(vvyF + 1) * (vvyF + 1) + (vvzF) * (vvzF)-1) +
+ // (mfbba) * ((vvxF) * (vvxF)+(vvyF) * (vvyF)+(vvzF + 1) * (vvzF + 1)-1) +
+ // (mfaab) * ((vvxF + 1) * (vvxF + 1) + (vvyF + 1) * (vvyF + 1) + (vvzF) * (vvzF)-1) +
+ // (mfcab) * ((vvxF - 1) * (vvxF - 1) + (vvyF + 1) * (vvyF + 1) + (vvzF) * (vvzF)-1) +
+ // (mfaba) * ((vvxF + 1) * (vvxF + 1) + (vvyF) * (vvyF)+(vvzF + 1) * (vvzF + 1)-1) +
+ // (mfcba) * ((vvxF - 1) * (vvxF - 1) + (vvyF) * (vvyF)+(vvzF + 1) * (vvzF + 1)-1) +
+ // (mfbaa) * ((vvxF) * (vvxF)+(vvyF + 1) * (vvyF + 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfbca) * ((vvxF) * (vvxF)+(vvyF - 1) * (vvyF - 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfaaa) * ((vvxF + 1) * (vvxF + 1) + (vvyF + 1) * (vvyF + 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfcaa) * ((vvxF - 1) * (vvxF - 1) + (vvyF + 1) * (vvyF + 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfaca) * ((vvxF + 1) * (vvxF + 1) + (vvyF - 1) * (vvyF - 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfcca) * ((vvxF - 1) * (vvxF - 1) + (vvyF - 1) * (vvyF - 1) + (vvzF + 1) * (vvzF + 1)-1) +
+ // (mfbbb) * ((vvxF) * (vvxF)+(vvyF) * (vvyF)+(vvzF) * (vvzF)-1);
+ //if (divAfterSource - divBeforeSource != 0 && phi[REST]>0.0001 && phi[REST]<0.999) {
+ // std::cout << phi[REST]<<" "<< divAfterSource << " " << divBeforeSource <<" "<< divAfterSource/ divBeforeSource << std::endl;
+ //}
+
+ //if (fabs(divAfterSource - divBeforeSource)/(fabs(divAfterSource) + fabs(divBeforeSource)+1e-10) > 1e-5) {
+ // LBMReal scaleDiv =0.95+(1-0.95)* (divBeforeSource) / (divBeforeSource - divAfterSource);
+
+ // forcingTerm[E] *=scaleDiv;
+ // forcingTerm[N] *=scaleDiv;
+ // forcingTerm[T] *=scaleDiv;
+ // forcingTerm[NE] *=scaleDiv;
+ // forcingTerm[NW] *=scaleDiv;
+ // forcingTerm[TE] *=scaleDiv;
+ // forcingTerm[TW] *=scaleDiv;
+ // forcingTerm[TN] *=scaleDiv;
+ // forcingTerm[TS] *=scaleDiv;
+ // forcingTerm[TNE] *=scaleDiv;
+ // forcingTerm[TNW] *=scaleDiv;
+ // forcingTerm[TSE] *=scaleDiv;
+ // forcingTerm[TSW] *=scaleDiv;
+ // forcingTerm[W] *=scaleDiv;
+ // forcingTerm[S] *=scaleDiv;
+ // forcingTerm[B] *=scaleDiv;
+ // forcingTerm[SW] *=scaleDiv;
+ // forcingTerm[SE] *=scaleDiv;
+ // forcingTerm[BW] *=scaleDiv;
+ // forcingTerm[BE] *=scaleDiv;
+ // forcingTerm[BS] *=scaleDiv;
+ // forcingTerm[BN] *=scaleDiv;
+ // forcingTerm[BSW] *=scaleDiv;
+ // forcingTerm[BSE] *=scaleDiv;
+ // forcingTerm[BNW] *=scaleDiv;
+ // forcingTerm[BNE] *=scaleDiv;
+ // forcingTerm[REST] *=scaleDiv;
+ //}
+ ////////
+
+
+ mfcbb += 3.0 * (0.5 * forcingTerm[E]) / rho; //-(3.0*p1 - rho)*WEIGTH[E ];
+ mfbcb += 3.0 * (0.5 * forcingTerm[N]) / rho; //-(3.0*p1 - rho)*WEIGTH[N ];
+ mfbbc += 3.0 * (0.5 * forcingTerm[T]) / rho; //-(3.0*p1 - rho)*WEIGTH[T ];
+ mfccb += 3.0 * (0.5 * forcingTerm[NE]) / rho; //-(3.0*p1 - rho)*WEIGTH[NE ];
+ mfacb += 3.0 * (0.5 * forcingTerm[NW]) / rho; //-(3.0*p1 - rho)*WEIGTH[NW ];
+ mfcbc += 3.0 * (0.5 * forcingTerm[TE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TE ];
+ mfabc += 3.0 * (0.5 * forcingTerm[TW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TW ];
+ mfbcc += 3.0 * (0.5 * forcingTerm[TN]) / rho; //-(3.0*p1 - rho)*WEIGTH[TN ];
+ mfbac += 3.0 * (0.5 * forcingTerm[TS]) / rho; //-(3.0*p1 - rho)*WEIGTH[TS ];
+ mfccc += 3.0 * (0.5 * forcingTerm[TNE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TNE];
+ mfacc += 3.0 * (0.5 * forcingTerm[TNW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TNW];
+ mfcac += 3.0 * (0.5 * forcingTerm[TSE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TSE];
+ mfaac += 3.0 * (0.5 * forcingTerm[TSW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TSW];
+ mfabb += 3.0 * (0.5 * forcingTerm[W]) / rho; //-(3.0*p1 - rho)*WEIGTH[W ];
+ mfbab += 3.0 * (0.5 * forcingTerm[S]) / rho; //-(3.0*p1 - rho)*WEIGTH[S ];
+ mfbba += 3.0 * (0.5 * forcingTerm[B]) / rho; //-(3.0*p1 - rho)*WEIGTH[B ];
+ mfaab += 3.0 * (0.5 * forcingTerm[SW]) / rho; //-(3.0*p1 - rho)*WEIGTH[SW ];
+ mfcab += 3.0 * (0.5 * forcingTerm[SE]) / rho; //-(3.0*p1 - rho)*WEIGTH[SE ];
+ mfaba += 3.0 * (0.5 * forcingTerm[BW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BW ];
+ mfcba += 3.0 * (0.5 * forcingTerm[BE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BE ];
+ mfbaa += 3.0 * (0.5 * forcingTerm[BS]) / rho; //-(3.0*p1 - rho)*WEIGTH[BS ];
+ mfbca += 3.0 * (0.5 * forcingTerm[BN]) / rho; //-(3.0*p1 - rho)*WEIGTH[BN ];
+ mfaaa += 3.0 * (0.5 * forcingTerm[BSW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BSW];
+ mfcaa += 3.0 * (0.5 * forcingTerm[BSE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BSE];
+ mfaca += 3.0 * (0.5 * forcingTerm[BNW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BNW];
+ mfcca += 3.0 * (0.5 * forcingTerm[BNE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BNE];
+ mfbbb += 3.0 * (0.5 * forcingTerm[REST]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
+
+ //--------------------------------------------------------
+
+
+
+
+
+ //////end classic source term
+
+
+
+
//forcing
///////////////////////////////////////////////////////////////////////////////////////////
if (withForcing)
@@ -418,13 +647,13 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
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();
+ //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
+ //vvx += forcingX1 * deltaT * 0.5; // X
+ //vvy += forcingX2 * deltaT * 0.5; // Y
+ //vvz += forcingX3 * deltaT * 0.5; // Z
}
LBMReal vx2;
@@ -697,12 +926,24 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
// Cumulants
////////////////////////////////////////////////////////////////////////////////////
LBMReal OxxPyyPzz = 1.; //omega2 or bulk viscosity
- LBMReal OxyyPxzz = 1.;//-s9;//2+s9;//
- LBMReal OxyyMxzz = 1.;//2+s9;//
+ // LBMReal OxyyPxzz = 1.;//-s9;//2+s9;//
+ // LBMReal OxyyMxzz = 1.;//2+s9;//
LBMReal O4 = 1.;
LBMReal O5 = 1.;
LBMReal O6 = 1.;
+
+
+ /////fourth order parameters; here only for test. Move out of loop!
+
+ LBMReal OxyyPxzz = 8.0 * (collFactorM - 2.0) * (OxxPyyPzz * (3.0 * collFactorM - 1.0) - 5.0 * collFactorM) / (8.0 * (5.0 - 2.0 * collFactorM) * collFactorM + OxxPyyPzz * (8.0 + collFactorM * (9.0 * collFactorM - 26.0)));
+ LBMReal OxyyMxzz = 8.0 * (collFactorM - 2.0) * (collFactorM + OxxPyyPzz * (3.0 * collFactorM - 7.0)) / (OxxPyyPzz * (56.0 - 42.0 * collFactorM + 9.0 * collFactorM * collFactorM) - 8.0 * collFactorM);
+ LBMReal Oxyz = 24.0 * (collFactorM - 2.0) * (4.0 * collFactorM * collFactorM + collFactorM * OxxPyyPzz * (18.0 - 13.0 * collFactorM) + OxxPyyPzz * OxxPyyPzz * (2.0 + collFactorM * (6.0 * collFactorM - 11.0))) / (16.0 * collFactorM * collFactorM * (collFactorM - 6.0) - 2.0 * collFactorM * OxxPyyPzz * (216.0 + 5.0 * collFactorM * (9.0 * collFactorM - 46.0)) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (3.0 * collFactorM - 10.0) * (15.0 * collFactorM - 28.0) - 48.0));
+ LBMReal A = (4.0 * collFactorM * collFactorM + 2.0 * collFactorM * OxxPyyPzz * (collFactorM - 6.0) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (10.0 - 3.0 * collFactorM) - 4.0)) / ((collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
+ //FIXME: warning C4459: declaration of 'B' hides global declaration (message : see declaration of 'D3Q27System::B' )
+ LBMReal BB = (4.0 * collFactorM * OxxPyyPzz * (9.0 * collFactorM - 16.0) - 4.0 * collFactorM * collFactorM - 2.0 * OxxPyyPzz * OxxPyyPzz * (2.0 + 9.0 * collFactorM * (collFactorM - 2.0))) / (3.0 * (collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
+
+
//Cum 4.
//LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
//LBMReal CUMbcb = mfbcb - ((mfaca + c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015
@@ -744,19 +985,24 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
//applying phase field gradients first part:
// mxxPyyPzz += c2o3 * rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz);
// 17.03.2021 attempt for statililization by assymptotically vanishing bias
- LBMReal correctionScaling =
- rhoToPhi / rho; // +0.5;// (vx2 + vy2 + vz2) * 100;// +0.5;//(vx2 + vy2 + vz2)*1000;
- mxxPyyPzz += (1.0 / 6.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) *
- correctionScaling; // As in Hesam's code
- mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy) * correctionScaling;
- mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
- mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
- mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
- mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
-
- LBMReal dxux = 0.0;// -c1o2 * collFactorM * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz);
- LBMReal dyuy = 0.0;// dxux + collFactorM * c3o2 * mxxMyy;
- LBMReal dzuz = 0.0;// dxux + collFactorM * c3o2 * mxxMzz;
+ //LBMReal correctionScaling =
+ // rhoToPhi / rho; // +0.5;// (vx2 + vy2 + vz2) * 100;// +0.5;//(vx2 + vy2 + vz2)*1000;
+ //mxxPyyPzz += (1.0 / 6.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) *
+ // correctionScaling; // As in Hesam's code
+ //mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy) * correctionScaling;
+ //mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
+ //mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
+ //mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
+ //mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
+
+ LBMReal dxux = -c1o2 * collFactorM * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz);
+ LBMReal dyuy = dxux + collFactorM * c3o2 * mxxMyy;
+ LBMReal dzuz = dxux + collFactorM * c3o2 * mxxMzz;
+
+ LBMReal Dxy = -three * collFactorM * mfbba;
+ LBMReal Dxz = -three * collFactorM * mfbab;
+ LBMReal Dyz = -three * collFactorM * mfabb;
+
//relax
mxxPyyPzz += OxxPyyPzz * (/*mfaaa*/ - mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
@@ -769,16 +1015,16 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
//applying phase field gradients second part:
//mxxPyyPzz += c2o3 * rhoToPhi * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz);
- mxxPyyPzz += (1.0 / 6.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) *
- correctionScaling; // As in Hesam's code
- mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy) * correctionScaling;
- mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
- mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
- mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
- mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
+ //mxxPyyPzz += (1.0 / 6.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) *
+ // correctionScaling; // As in Hesam's code
+ //mxxMyy += c1o3 * (dX1_phi * vvx - dX2_phi * vvy) * correctionScaling;
+ //mxxMzz += c1o3 * (dX1_phi * vvx - dX3_phi * vvz) * correctionScaling;
+ //mfabb += c1o6 * (dX2_phi * vvz + dX3_phi * vvy) * correctionScaling;
+ //mfbab += c1o6 * (dX1_phi * vvz + dX3_phi * vvx) * correctionScaling;
+ //mfbba += c1o6 * (dX1_phi * vvy + dX2_phi * vvx) * correctionScaling;
////updated pressure
- mfaaa += (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling;
+ //mfaaa += (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling;
mxxPyyPzz += mfaaa; // 12.03.21 shifted by mfaaa
// linear combinations back
@@ -822,13 +1068,19 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
mfbac = (-mxyyMxzz + mxyyPxzz) * c1o2;
//4.
- CUMacc += O4 * (-CUMacc);
- CUMcac += O4 * (-CUMcac);
- CUMcca += O4 * (-CUMcca);
-
- CUMbbc += O4 * (-CUMbbc);
- CUMbcb += O4 * (-CUMbcb);
- CUMcbb += O4 * (-CUMcbb);
+ //CUMacc += O4 * (-CUMacc);
+ //CUMcac += O4 * (-CUMcac);
+ //CUMcca += O4 * (-CUMcca);
+
+ //CUMbbc += O4 * (-CUMbbc);
+ //CUMbcb += O4 * (-CUMbcb);
+ //CUMcbb += O4 * (-CUMcbb);
+ CUMacc = -O4 * (one / collFactorM - c1o2) * (dyuy + dzuz) * c2o3 * A + (one - O4) * (CUMacc);
+ CUMcac = -O4 * (one / collFactorM - c1o2) * (dxux + dzuz) * c2o3 * A + (one - O4) * (CUMcac);
+ CUMcca = -O4 * (one / collFactorM - c1o2) * (dyuy + dxux) * c2o3 * A + (one - O4) * (CUMcca);
+ CUMbbc = -O4 * (one / collFactorM - c1o2) * Dxy * c1o3 * BB + (one - O4) * (CUMbbc);
+ CUMbcb = -O4 * (one / collFactorM - c1o2) * Dxz * c1o3 * BB + (one - O4) * (CUMbcb);
+ CUMcbb = -O4 * (one / collFactorM - c1o2) * Dyz * c1o3 * BB + (one - O4) * (CUMcbb);
//5.
CUMbcc += O5 * (-CUMbcc);
@@ -1090,60 +1342,92 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
mfbcc = m1;
mfccc = m2;
+ /////classical source term 8.4.2021
+
+ mfcbb += 3.0 * (0.5 * forcingTerm[E]) / rho; //-(3.0*p1 - rho)*WEIGTH[E ];
+ mfbcb += 3.0 * (0.5 * forcingTerm[N]) / rho; //-(3.0*p1 - rho)*WEIGTH[N ];
+ mfbbc += 3.0 * (0.5 * forcingTerm[T]) / rho; //-(3.0*p1 - rho)*WEIGTH[T ];
+ mfccb += 3.0 * (0.5 * forcingTerm[NE]) / rho; //-(3.0*p1 - rho)*WEIGTH[NE ];
+ mfacb += 3.0 * (0.5 * forcingTerm[NW]) / rho; //-(3.0*p1 - rho)*WEIGTH[NW ];
+ mfcbc += 3.0 * (0.5 * forcingTerm[TE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TE ];
+ mfabc += 3.0 * (0.5 * forcingTerm[TW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TW ];
+ mfbcc += 3.0 * (0.5 * forcingTerm[TN]) / rho; //-(3.0*p1 - rho)*WEIGTH[TN ];
+ mfbac += 3.0 * (0.5 * forcingTerm[TS]) / rho; //-(3.0*p1 - rho)*WEIGTH[TS ];
+ mfccc += 3.0 * (0.5 * forcingTerm[TNE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TNE];
+ mfacc += 3.0 * (0.5 * forcingTerm[TNW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TNW];
+ mfcac += 3.0 * (0.5 * forcingTerm[TSE]) / rho; //-(3.0*p1 - rho)*WEIGTH[TSE];
+ mfaac += 3.0 * (0.5 * forcingTerm[TSW]) / rho; //-(3.0*p1 - rho)*WEIGTH[TSW];
+ mfabb += 3.0 * (0.5 * forcingTerm[W]) / rho; //-(3.0*p1 - rho)*WEIGTH[W ];
+ mfbab += 3.0 * (0.5 * forcingTerm[S]) / rho; //-(3.0*p1 - rho)*WEIGTH[S ];
+ mfbba += 3.0 * (0.5 * forcingTerm[B]) / rho; //-(3.0*p1 - rho)*WEIGTH[B ];
+ mfaab += 3.0 * (0.5 * forcingTerm[SW]) / rho; //-(3.0*p1 - rho)*WEIGTH[SW ];
+ mfcab += 3.0 * (0.5 * forcingTerm[SE]) / rho; //-(3.0*p1 - rho)*WEIGTH[SE ];
+ mfaba += 3.0 * (0.5 * forcingTerm[BW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BW ];
+ mfcba += 3.0 * (0.5 * forcingTerm[BE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BE ];
+ mfbaa += 3.0 * (0.5 * forcingTerm[BS]) / rho; //-(3.0*p1 - rho)*WEIGTH[BS ];
+ mfbca += 3.0 * (0.5 * forcingTerm[BN]) / rho; //-(3.0*p1 - rho)*WEIGTH[BN ];
+ mfaaa += 3.0 * (0.5 * forcingTerm[BSW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BSW];
+ mfcaa += 3.0 * (0.5 * forcingTerm[BSE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BSE];
+ mfaca += 3.0 * (0.5 * forcingTerm[BNW]) / rho; //-(3.0*p1 - rho)*WEIGTH[BNW];
+ mfcca += 3.0 * (0.5 * forcingTerm[BNE]) / rho; //-(3.0*p1 - rho)*WEIGTH[BNE];
+ mfbbb += 3.0 * (0.5 * forcingTerm[REST]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
+
+
+
//////////////////////////////////////////////////////////////////////////
//proof correctness
//////////////////////////////////////////////////////////////////////////
-#ifdef PROOF_CORRECTNESS
- LBMReal rho_post = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca)
- + (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) + (mfbaa + mfbac + mfbca + mfbcc)
- + (mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
- //LBMReal dif = fabs(drho - rho_post);
- LBMReal dif = drho + (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling - rho_post;
-#ifdef SINGLEPRECISION
- if (dif > 10.0E-7 || dif < -10.0E-7)
-#else
- if (dif > 10.0E-15 || dif < -10.0E-15)
-#endif
- {
- UB_THROW(UbException(UB_EXARGS, "drho=" + UbSystem::toString(drho) + ", rho_post=" + UbSystem::toString(rho_post)
- + " dif=" + UbSystem::toString(dif)
- + " drho is not correct for node " + UbSystem::toString(x1) + "," + UbSystem::toString(x2) + "," + UbSystem::toString(x3)));
- //UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): drho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3));
- //exit(EXIT_FAILURE);
- }
-#endif
+//#ifdef PROOF_CORRECTNESS
+// LBMReal rho_post = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca)
+// + (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) + (mfbaa + mfbac + mfbca + mfbcc)
+// + (mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
+// //LBMReal dif = fabs(drho - rho_post);
+// LBMReal dif = drho + (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling - rho_post;
+//#ifdef SINGLEPRECISION
+// if (dif > 10.0E-7 || dif < -10.0E-7)
+//#else
+// if (dif > 10.0E-15 || dif < -10.0E-15)
+//#endif
+// {
+// UB_THROW(UbException(UB_EXARGS, "drho=" + UbSystem::toString(drho) + ", rho_post=" + UbSystem::toString(rho_post)
+// + " dif=" + UbSystem::toString(dif)
+// + " drho is not correct for node " + UbSystem::toString(x1) + "," + UbSystem::toString(x2) + "," + UbSystem::toString(x3)));
+// //UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): drho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3));
+// //exit(EXIT_FAILURE);
+// }
+//#endif
//////////////////////////////////////////////////////////////////////////
//write distribution
//////////////////////////////////////////////////////////////////////////
- (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) = mfabb * rho;
- (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) = mfbab * rho;
- (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3) = mfbba * rho;
- (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3) = mfaab * rho;
- (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3) = mfcab * rho;
- (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3) = mfaba * rho;
- (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3) = mfcba * rho;
- (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3) = mfbaa * rho;
- (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3) = mfbca * rho;
- (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3) = mfaaa * rho;
- (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3) = mfcaa * rho;
- (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3) = mfaca * rho;
- (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfcca * rho;
-
- (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3) = mfcbb * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3) = mfbcb * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p) = mfbbc * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3) = mfccb * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3) = mfacb * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p) = mfcbc * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p) = mfabc * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbcc * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p) = mfbac * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfccc * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfacc * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac * rho;
- (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac * rho;
-
- (*this->zeroDistributionsF)(x1, x2, x3) = mfbbb * rho;
+ (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) = mfabb * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) = mfbab * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3) = mfbba * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3) = mfaab * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3) = mfcab * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3) = mfaba * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3) = mfcba * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3) = mfbaa * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3) = mfbca * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3) = mfaaa * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3) = mfcaa * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3) = mfaca * rho * c1o3;
+ (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfcca * rho * c1o3;
+
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3) = mfcbb * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3) = mfbcb * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p) = mfbbc * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3) = mfccb * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3) = mfacb * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p) = mfcbc * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p) = mfabc * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbcc * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p) = mfbac * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfccc * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfacc * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac * rho * c1o3;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac * rho * c1o3;
+
+ (*this->zeroDistributionsF)(x1, x2, x3) = mfbbb * rho * c1o3;
//////////////////////////////////////////////////////////////////////////
////!Incompressible Kernal
@@ -1977,7 +2261,7 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
////////////////////////////////////////////
/////CUMULANT PHASE-FIELD
LBMReal omegaD =1.0/( 3.0 * mob + 0.5);
-
+ {
mfcbb = (*this->localDistributionsH1)(D3Q27System::ET_E, x1, x2, x3);
mfbcb = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
mfbbc = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
@@ -2134,6 +2418,9 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
LBMReal Mccb = mfccb - mfaab * c1o9;
// collision of 1st order moments
+ // LBMReal ccx, ccy, ccz;
+
+
cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
@@ -2281,6 +2568,323 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
(*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac;
(*this->zeroDistributionsH1)(x1,x2,x3) = mfbbb;
+ }
+
+ ////Phasefield 2:
+
+ {
+
+ normX1 *= -1;
+ normX2 *= -1;
+ normX3 *= -1;
+
+ mfcbb = (*this->localDistributionsH2)(D3Q27System::ET_E, x1, x2, x3);
+ mfbcb = (*this->localDistributionsH2)(D3Q27System::ET_N, x1, x2, x3);
+ mfbbc = (*this->localDistributionsH2)(D3Q27System::ET_T, x1, x2, x3);
+ mfccb = (*this->localDistributionsH2)(D3Q27System::ET_NE, x1, x2, x3);
+ mfacb = (*this->localDistributionsH2)(D3Q27System::ET_NW, x1p, x2, x3);
+ mfcbc = (*this->localDistributionsH2)(D3Q27System::ET_TE, x1, x2, x3);
+ mfabc = (*this->localDistributionsH2)(D3Q27System::ET_TW, x1p, x2, x3);
+ mfbcc = (*this->localDistributionsH2)(D3Q27System::ET_TN, x1, x2, x3);
+ mfbac = (*this->localDistributionsH2)(D3Q27System::ET_TS, x1, x2p, x3);
+ mfccc = (*this->localDistributionsH2)(D3Q27System::ET_TNE, x1, x2, x3);
+ mfacc = (*this->localDistributionsH2)(D3Q27System::ET_TNW, x1p, x2, x3);
+ mfcac = (*this->localDistributionsH2)(D3Q27System::ET_TSE, x1, x2p, x3);
+ mfaac = (*this->localDistributionsH2)(D3Q27System::ET_TSW, x1p, x2p, x3);
+ mfabb = (*this->nonLocalDistributionsH2)(D3Q27System::ET_W, x1p, x2, x3);
+ mfbab = (*this->nonLocalDistributionsH2)(D3Q27System::ET_S, x1, x2p, x3);
+ mfbba = (*this->nonLocalDistributionsH2)(D3Q27System::ET_B, x1, x2, x3p);
+ mfaab = (*this->nonLocalDistributionsH2)(D3Q27System::ET_SW, x1p, x2p, x3);
+ mfcab = (*this->nonLocalDistributionsH2)(D3Q27System::ET_SE, x1, x2p, x3);
+ mfaba = (*this->nonLocalDistributionsH2)(D3Q27System::ET_BW, x1p, x2, x3p);
+ mfcba = (*this->nonLocalDistributionsH2)(D3Q27System::ET_BE, x1, x2, x3p);
+ mfbaa = (*this->nonLocalDistributionsH2)(D3Q27System::ET_BS, x1, x2p, x3p);
+ mfbca = (*this->nonLocalDistributionsH2)(D3Q27System::ET_BN, x1, x2, x3p);
+ mfaaa = (*this->nonLocalDistributionsH2)(D3Q27System::ET_BSW, x1p, x2p, x3p);
+ mfcaa = (*this->nonLocalDistributionsH2)(D3Q27System::ET_BSE, x1, x2p, x3p);
+ mfaca = (*this->nonLocalDistributionsH2)(D3Q27System::ET_BNW, x1p, x2, x3p);
+ mfcca = (*this->nonLocalDistributionsH2)(D3Q27System::ET_BNE, x1, x2, x3p);
+ mfbbb = (*this->zeroDistributionsH2)(x1, x2, x3);
+
+
+ ////////////////////////////////////////////////////////////////////////////////////
+//! - Calculate density and velocity using pyramid summation for low round-off errors as in Eq. (J1)-(J3) \ref
+//! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
+//!
+////////////////////////////////////////////////////////////////////////////////////
+// fluid component
+ //LBMReal drhoFluid =
+ // ((((fccc + faaa) + (faca + fcac)) + ((facc + fcaa) + (faac + fcca))) +
+ // (((fbac + fbca) + (fbaa + fbcc)) + ((fabc + fcba) + (faba + fcbc)) + ((facb + fcab) + (faab + fccb))) +
+ // ((fabb + fcbb) + (fbab + fbcb) + (fbba + fbbc))) + fbbb;
+
+ //LBMReal rhoFluid = c1 + drhoFluid;
+ //LBMReal OOrhoFluid = c1 / rhoFluid;
+
+
+ //LBMReal vvx =
+ // ((((fccc - faaa) + (fcac - faca)) + ((fcaa - facc) + (fcca - faac))) +
+ // (((fcba - fabc) + (fcbc - faba)) + ((fcab - facb) + (fccb - faab))) +
+ // (fcbb - fabb)) * OOrhoFluid;
+ //LBMReal vvy =
+ // ((((fccc - faaa) + (faca - fcac)) + ((facc - fcaa) + (fcca - faac))) +
+ // (((fbca - fbac) + (fbcc - fbaa)) + ((facb - fcab) + (fccb - faab))) +
+ // (fbcb - fbab)) * OOrhoFluid;
+ //LBMReal vvz =
+ // ((((fccc - faaa) + (fcac - faca)) + ((facc - fcaa) + (faac - fcca))) +
+ // (((fbac - fbca) + (fbcc - fbaa)) + ((fabc - fcba) + (fcbc - faba))) +
+ // (fbbc - fbba)) * OOrhoFluid;
+
+ // LBMReal vvx = ux;
+ // LBMReal vvy = uy;
+ // LBMReal vvz = uz;
+ ////////////////////////////////////////////////////////////////////////////////////
+ // second component
+ LBMReal concentration =
+ ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) +
+ (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
+ ((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Add half of the acceleration (body force) to the velocity as in Eq. (42) \ref
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
+ //!
+ // LBMReal fx = forces[0];
+ // LBMReal fy = forces[1];
+ // LBMReal fz = -concentration * forces[2];
+ // vvx += fx * c1o2;
+ // vvy += fy * c1o2;
+ // vvz += fz * c1o2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ LBMReal oneMinusRho = c1 - concentration;
+
+ LBMReal cx =
+ ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+ (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
+ (mfcbb - mfabb));
+ LBMReal cy =
+ ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+ (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
+ (mfbcb - mfbab));
+ LBMReal cz =
+ ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+ (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
+ (mfbbc - mfbba));
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // calculate the square of velocities for this lattice node
+ LBMReal cx2 = cx * cx;
+ LBMReal cy2 = cy * cy;
+ LBMReal cz2 = cz * cz;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Chimera transform from well conditioned distributions to central moments as defined in Appendix J in \ref
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
+ //! see also Eq. (6)-(14) in \ref
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ forwardInverseChimeraWithKincompressible(mfaaa, mfaab, mfaac, cz, cz2, c36, c1o36, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfaba, mfabb, mfabc, cz, cz2, c9, c1o9, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfaca, mfacb, mfacc, cz, cz2, c36, c1o36, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfbaa, mfbab, mfbac, cz, cz2, c9, c1o9, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfbba, mfbbb, mfbbc, cz, cz2, c9o4, c4o9, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfbca, mfbcb, mfbcc, cz, cz2, c9, c1o9, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfcaa, mfcab, mfcac, cz, cz2, c36, c1o36, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfcba, mfcbb, mfcbc, cz, cz2, c9, c1o9, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfcca, mfccb, mfccc, cz, cz2, c36, c1o36, oneMinusRho);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ forwardInverseChimeraWithKincompressible(mfaaa, mfaba, mfaca, cy, cy2, c6, c1o6, oneMinusRho);
+ forwardChimera(mfaab, mfabb, mfacb, cy, cy2);
+ forwardInverseChimeraWithKincompressible(mfaac, mfabc, mfacc, cy, cy2, c18, c1o18, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfbaa, mfbba, mfbca, cy, cy2, c3o2, c2o3, oneMinusRho);
+ forwardChimera(mfbab, mfbbb, mfbcb, cy, cy2);
+ forwardInverseChimeraWithKincompressible(mfbac, mfbbc, mfbcc, cy, cy2, c9o2, c2o9, oneMinusRho);
+ forwardInverseChimeraWithKincompressible(mfcaa, mfcba, mfcca, cy, cy2, c6, c1o6, oneMinusRho);
+ forwardChimera(mfcab, mfcbb, mfccb, cy, cy2);
+ forwardInverseChimeraWithKincompressible(mfcac, mfcbc, mfccc, cy, cy2, c18, c1o18, oneMinusRho);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ forwardInverseChimeraWithKincompressible(mfaaa, mfbaa, mfcaa, cx, cx2, c1, c1, oneMinusRho);
+ forwardChimera(mfaba, mfbba, mfcba, cx, cx2);
+ forwardInverseChimeraWithKincompressible(mfaca, mfbca, mfcca, cx, cx2, c3, c1o3, oneMinusRho);
+ forwardChimera(mfaab, mfbab, mfcab, cx, cx2);
+ forwardChimera(mfabb, mfbbb, mfcbb, cx, cx2);
+ forwardChimera(mfacb, mfbcb, mfccb, cx, cx2);
+ forwardInverseChimeraWithKincompressible(mfaac, mfbac, mfcac, cx, cx2, c3, c1o3, oneMinusRho);
+ forwardChimera(mfabc, mfbbc, mfcbc, cx, cx2);
+ forwardInverseChimeraWithKincompressible(mfacc, mfbcc, mfccc, cx, cx2, c3, c1o9, oneMinusRho);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - experimental Cumulant ... to be published ... hopefully
+ //!
+
+ // linearized orthogonalization of 3rd order central moments
+ LBMReal Mabc = mfabc - mfaba * c1o3;
+ LBMReal Mbca = mfbca - mfbaa * c1o3;
+ LBMReal Macb = mfacb - mfaab * c1o3;
+ LBMReal Mcba = mfcba - mfaba * c1o3;
+ LBMReal Mcab = mfcab - mfaab * c1o3;
+ LBMReal Mbac = mfbac - mfbaa * c1o3;
+ // linearized orthogonalization of 5th order central moments
+ LBMReal Mcbc = mfcbc - mfaba * c1o9;
+ LBMReal Mbcc = mfbcc - mfbaa * c1o9;
+ LBMReal Mccb = mfccb - mfaab * c1o9;
+
+ // collision of 1st order moments
+ cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
+ normX1 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+ cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
+ normX2 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+ cz = cz * (c1 - omegaD) + omegaD * vvz * concentration +
+ normX3 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+
+ //mhx = (ux * phi[REST] + normX1 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhx;
+ //mhy = (uy * phi[REST] + normX2 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhy;
+//mhz = (uz * phi[REST] + normX3 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhz;
+
+
+ cx2 = cx * cx;
+ cy2 = cy * cy;
+ cz2 = cz * cz;
+
+ // equilibration of 2nd order moments
+ mfbba = zeroReal;
+ mfbab = zeroReal;
+ mfabb = zeroReal;
+
+ mfcaa = c1o3 * concentration;
+ mfaca = c1o3 * concentration;
+ mfaac = c1o3 * concentration;
+
+
+ //LBMReal omega2 = 1.0f;// omegaD;
+ //mfbba *= (c1 - omega2);
+ //mfbab *= (c1 - omega2);
+ //mfabb *= (c1 - omega2);
+
+ //mfcaa = mfcaa*(c1 - omega2) + omega2*c1o3 * concentration;
+ //mfaca = mfaca*(c1 - omega2) + omega2*c1o3 * concentration;
+ //mfaac = mfaac*(c1 - omega2) + omega2*c1o3 * concentration;
+
+ // equilibration of 3rd order moments
+ Mabc = zeroReal;
+ Mbca = zeroReal;
+ Macb = zeroReal;
+ Mcba = zeroReal;
+ Mcab = zeroReal;
+ Mbac = zeroReal;
+ mfbbb = zeroReal;
+
+ // from linearized orthogonalization 3rd order central moments to central moments
+ mfabc = Mabc + mfaba * c1o3;
+ mfbca = Mbca + mfbaa * c1o3;
+ mfacb = Macb + mfaab * c1o3;
+ mfcba = Mcba + mfaba * c1o3;
+ mfcab = Mcab + mfaab * c1o3;
+ mfbac = Mbac + mfbaa * c1o3;
+
+ // equilibration of 4th order moments
+ mfacc = c1o9 * concentration;
+ mfcac = c1o9 * concentration;
+ mfcca = c1o9 * concentration;
+
+ mfcbb = zeroReal;
+ mfbcb = zeroReal;
+ mfbbc = zeroReal;
+
+ // equilibration of 5th order moments
+ Mcbc = zeroReal;
+ Mbcc = zeroReal;
+ Mccb = zeroReal;
+
+ // from linearized orthogonalization 5th order central moments to central moments
+ mfcbc = Mcbc + mfaba * c1o9;
+ mfbcc = Mbcc + mfbaa * c1o9;
+ mfccb = Mccb + mfaab * c1o9;
+
+ // equilibration of 6th order moment
+ mfccc = c1o27 * concentration;
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Chimera transform from central moments to well conditioned distributions as defined in Appendix J in
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
+ //! see also Eq. (88)-(96) in
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ backwardInverseChimeraWithKincompressible(mfaaa, mfbaa, mfcaa, cx, cx2, c1, c1, oneMinusRho);
+ backwardChimera(mfaba, mfbba, mfcba, cx, cx2);
+ backwardInverseChimeraWithKincompressible(mfaca, mfbca, mfcca, cx, cx2, c3, c1o3, oneMinusRho);
+ backwardChimera(mfaab, mfbab, mfcab, cx, cx2);
+ backwardChimera(mfabb, mfbbb, mfcbb, cx, cx2);
+ backwardChimera(mfacb, mfbcb, mfccb, cx, cx2);
+ backwardInverseChimeraWithKincompressible(mfaac, mfbac, mfcac, cx, cx2, c3, c1o3, oneMinusRho);
+ backwardChimera(mfabc, mfbbc, mfcbc, cx, cx2);
+ backwardInverseChimeraWithKincompressible(mfacc, mfbcc, mfccc, cx, cx2, c9, c1o9, oneMinusRho);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ backwardInverseChimeraWithKincompressible(mfaaa, mfaba, mfaca, cy, cy2, c6, c1o6, oneMinusRho);
+ backwardChimera(mfaab, mfabb, mfacb, cy, cy2);
+ backwardInverseChimeraWithKincompressible(mfaac, mfabc, mfacc, cy, cy2, c18, c1o18, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfbaa, mfbba, mfbca, cy, cy2, c3o2, c2o3, oneMinusRho);
+ backwardChimera(mfbab, mfbbb, mfbcb, cy, cy2);
+ backwardInverseChimeraWithKincompressible(mfbac, mfbbc, mfbcc, cy, cy2, c9o2, c2o9, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfcaa, mfcba, mfcca, cy, cy2, c6, c1o6, oneMinusRho);
+ backwardChimera(mfcab, mfcbb, mfccb, cy, cy2);
+ backwardInverseChimeraWithKincompressible(mfcac, mfcbc, mfccc, cy, cy2, c18, c1o18, oneMinusRho);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ backwardInverseChimeraWithKincompressible(mfaaa, mfaab, mfaac, cz, cz2, c36, c1o36, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfaba, mfabb, mfabc, cz, cz2, c9, c1o9, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfaca, mfacb, mfacc, cz, cz2, c36, c1o36, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfbaa, mfbab, mfbac, cz, cz2, c9, c1o9, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfbba, mfbbb, mfbbc, cz, cz2, c9o4, c4o9, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfbca, mfbcb, mfbcc, cz, cz2, c9, c1o9, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfcaa, mfcab, mfcac, cz, cz2, c36, c1o36, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfcba, mfcbb, mfcbc, cz, cz2, c9, c1o9, oneMinusRho);
+ backwardInverseChimeraWithKincompressible(mfcca, mfccb, mfccc, cz, cz2, c36, c1o36, oneMinusRho);
+
+
+
+ (*this->localDistributionsH2)(D3Q27System::ET_E, x1, x2, x3) = mfabb;
+ (*this->localDistributionsH2)(D3Q27System::ET_N, x1, x2, x3) = mfbab;
+ (*this->localDistributionsH2)(D3Q27System::ET_T, x1, x2, x3) = mfbba;
+ (*this->localDistributionsH2)(D3Q27System::ET_NE, x1, x2, x3) = mfaab;
+ (*this->localDistributionsH2)(D3Q27System::ET_NW, x1p, x2, x3) = mfcab;
+ (*this->localDistributionsH2)(D3Q27System::ET_TE, x1, x2, x3) = mfaba;
+ (*this->localDistributionsH2)(D3Q27System::ET_TW, x1p, x2, x3) = mfcba;
+ (*this->localDistributionsH2)(D3Q27System::ET_TN, x1, x2, x3) = mfbaa;
+ (*this->localDistributionsH2)(D3Q27System::ET_TS, x1, x2p, x3) = mfbca;
+ (*this->localDistributionsH2)(D3Q27System::ET_TNE, x1, x2, x3) = mfaaa;
+ (*this->localDistributionsH2)(D3Q27System::ET_TNW, x1p, x2, x3) = mfcaa;
+ (*this->localDistributionsH2)(D3Q27System::ET_TSE, x1, x2p, x3) = mfaca;
+ (*this->localDistributionsH2)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfcca;
+
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_W, x1p, x2, x3) = mfcbb;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_S, x1, x2p, x3) = mfbcb;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_B, x1, x2, x3p) = mfbbc;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_SW, x1p, x2p, x3) = mfccb;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_SE, x1, x2p, x3) = mfacb;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_BW, x1p, x2, x3p) = mfcbc;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_BE, x1, x2, x3p) = mfabc;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbcc;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_BN, x1, x2, x3p) = mfbac;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfccc;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfacc;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac;
+ (*this->nonLocalDistributionsH2)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac;
+
+ (*this->zeroDistributionsH2)(x1, x2, x3) = mfbbb;
+
+ }
+
+
/////!CUMULANT PHASE-FIELD
@@ -2371,8 +2975,8 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
}
}
}
- }
- dataSet->setPhaseField(divU);
+
+ // dataSet->setPhaseField(divU);
}
}
//////////////////////////////////////////////////////////////////////////
@@ -2416,6 +3020,49 @@ LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX3_phi()
//return 3.0 * sum;
}
+LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX1_phi2()
+{
+ using namespace D3Q27System;
+ return 3.0 * ((WEIGTH[TNE] * (((phi2[TNE] - phi2[BSW]) + (phi2[BSE] - phi2[TNW])) + ((phi2[TSE] - phi2[BNW]) + (phi2[BNE] - phi2[TSW])))
+ + WEIGTH[NE] * (((phi2[TE] - phi2[BW]) + (phi2[BE] - phi2[TW])) + ((phi2[SE] - phi2[NW]) + (phi2[NE] - phi2[SW])))) +
+ +WEIGTH[N] * (phi2[E] - phi2[W]));
+ //LBMReal sum = 0.0;
+ //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+ // sum += WEIGTH[k] * DX1[k] * phi2[k];
+ //}
+ //return 3.0 * sum;
+}
+
+LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX2_phi2()
+{
+ using namespace D3Q27System;
+ return 3.0 * ((WEIGTH[TNE] * (((phi2[TNE] - phi2[BSW]) - (phi2[BSE] - phi2[TNW])) + ((phi2[BNE] - phi2[TSW]) - (phi2[TSE] - phi2[BNW])))
+ + WEIGTH[NE] * (((phi2[TN] - phi2[BS]) + (phi2[BN] - phi2[TS])) + ((phi2[NE] - phi2[SW]) - (phi2[SE] - phi2[NW])))) +
+ +WEIGTH[N] * (phi2[N] - phi2[S]));
+ //LBMReal sum = 0.0;
+ //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+ // sum += WEIGTH[k] * DX2[k] * phi2[k];
+ //}
+ //return 3.0 * sum;
+}
+
+LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX3_phi2()
+{
+ using namespace D3Q27System;
+ return 3.0 * ((WEIGTH[TNE] * (((phi2[TNE] - phi2[BSW]) - (phi2[BSE] - phi2[TNW])) + ((phi2[TSE] - phi2[BNW]) - (phi2[BNE] - phi2[TSW])))
+ + WEIGTH[NE] * (((phi2[TE] - phi2[BW]) - (phi2[BE] - phi2[TW])) + ((phi2[TS] - phi2[BN]) + (phi2[TN] - phi2[BS])))) +
+ +WEIGTH[N] * (phi2[T] - phi2[B]));
+ //LBMReal sum = 0.0;
+ //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+ // sum += WEIGTH[k] * DX3[k] * phi2[k];
+ //}
+ //return 3.0 * sum;
+}
+
+
+
+
+
LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::nabla2_phi()
{
using namespace D3Q27System;
@@ -2504,6 +3151,7 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::findNeighbors(CbArray3D<LBMReal,
phi[REST] = (*ph)(x1, x2, x3);
+
for (int k = FSTARTDIR; k <= FENDDIR; k++) {
if (!bcArray->isSolid(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k])) {
@@ -2514,8 +3162,30 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::findNeighbors(CbArray3D<LBMReal,
}
}
+void MultiphaseTwoPhaseFieldsCumulantLBMKernel::findNeighbors2(CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr ph, int x1, int x2,
+ int x3)
+{
+ using namespace D3Q27System;
+
+ SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
+
+ phi2[REST] = (*ph)(x1, x2, x3);
+
+
+ for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+
+ if (!bcArray->isSolid(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k])) {
+ phi2[k] = (*ph)(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k]);
+ }
+ else {
+ phi2[k] = 0.0;
+ }
+ }
+}
+
void MultiphaseTwoPhaseFieldsCumulantLBMKernel::swapDistributions()
{
LBMKernel::swapDistributions();
dataSet->getHdistributions()->swap();
+ dataSet->getH2distributions()->swap();
}
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.h
index bc7609a29..a65fe073f 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.h
@@ -80,18 +80,23 @@ protected:
LBMReal h2[D3Q27System::ENDF + 1];
LBMReal g [D3Q27System::ENDF+1];
LBMReal phi[D3Q27System::ENDF+1];
+ LBMReal phi2[D3Q27System::ENDF + 1];
LBMReal pr1[D3Q27System::ENDF+1];
LBMReal phi_cutoff[D3Q27System::ENDF+1];
LBMReal gradX1_phi();
LBMReal gradX2_phi();
LBMReal gradX3_phi();
+ LBMReal gradX1_phi2();
+ LBMReal gradX2_phi2();
+ LBMReal gradX3_phi2();
//LBMReal gradX1_pr1();
//LBMReal gradX2_pr1();
//LBMReal gradX3_pr1();
//LBMReal dirgradC_phi(int n, int k);
void computePhasefield();
void findNeighbors(CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr ph /*Phase-Field*/, int x1, int x2, int x3);
+ void findNeighbors2(CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr ph, int x1, int x2, int x3);
//void findNeighbors(CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr ph /*Phase-Field*/, CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr pf /*Pressure-Field*/, int x1, int x2, int x3);
//void pressureFiltering(CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr pf /*Pressure-Field*/, CbArray3D<LBMReal,IndexerX3X2X1>::CbArray3DPtr pf_filtered /*Pressure-Field*/);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
index 4213e2b79..124e342d6 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
@@ -224,7 +224,7 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
feq[dir] = rho*WEIGTH[dir]*(1 + 3*velProd + 4.5*velSq1 - 1.5*(vx1Sq+vx2Sq+vx3Sq));
//geq[dir] = p1*WEIGTH1[dir] + gamma;
- geq[dir] = p1*WEIGTH[dir]/(rho*UbMath::c1o3) + gamma;
+ geq[dir] = p1*WEIGTH[dir]/(rho*UbMath::c1o3) + gamma*rho;
}
@@ -291,6 +291,35 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
distributionsH->setDistributionInv(f, ix1, ix2, ix3);
if (distributionsH2) {
+
+ f[E] = (1.-phi) * feq[E] / rho;
+ f[W] = (1.-phi) * feq[W] / rho;
+ f[N] = (1.-phi) * feq[N] / rho;
+ f[S] = (1.-phi) * feq[S] / rho;
+ f[T] = (1.-phi) * feq[T] / rho;
+ f[B] = (1.-phi) * feq[B] / rho;
+ f[NE] = (1.-phi) * feq[NE] / rho;
+ f[SW] = (1.-phi) * feq[SW] / rho;
+ f[SE] = (1.-phi) * feq[SE] / rho;
+ f[NW] = (1.-phi) * feq[NW] / rho;
+ f[TE] = (1.-phi) * feq[TE] / rho;
+ f[BW] = (1.-phi) * feq[BW] / rho;
+ f[BE] = (1.-phi) * feq[BE] / rho;
+ f[TW] = (1.-phi) * feq[TW] / rho;
+ f[TN] = (1.-phi) * feq[TN] / rho;
+ f[BS] = (1.-phi) * feq[BS] / rho;
+ f[BN] = (1.-phi) * feq[BN] / rho;
+ f[TS] = (1.-phi) * feq[TS] / rho;
+ f[TNE] = (1.-phi) * feq[TNE] / rho;
+ f[TNW] = (1.-phi) * feq[TNW] / rho;
+ f[TSE] = (1.-phi) * feq[TSE] / rho;
+ f[TSW] = (1.-phi) * feq[TSW] / rho;
+ f[BNE] = (1.-phi) * feq[BNE] / rho;
+ f[BNW] = (1.-phi) * feq[BNW] / rho;
+ f[BSE] = (1.-phi) * feq[BSE] / rho;
+ f[BSW] = (1.-phi) * feq[BSW] / rho;
+ f[REST] = (1.-phi) * feq[REST] / rho;
+
distributionsH2->setDistribution(f, ix1, ix2, ix3);
distributionsH2->setDistributionInv(f, ix1, ix2, ix3);
}
--
GitLab