diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
index 09c43e97b9b9f4ba46fb741b5f2cc2b1299d109c..f219a432a0535e4045cda7ef224485b97efc4282 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
@@ -85,49 +85,18 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
using namespace D3Q27System;
using namespace UbMath;
- // initializing of forcing stuff
- /*if (withForcing)
- {
- muForcingX1.DefineVar("x1",&muX1); muForcingX1.DefineVar("x2",&muX2); muForcingX1.DefineVar("x3",&muX3);
- muForcingX2.DefineVar("x1",&muX1); muForcingX2.DefineVar("x2",&muX2); muForcingX2.DefineVar("x3",&muX3);
- muForcingX3.DefineVar("x1",&muX1); muForcingX3.DefineVar("x2",&muX2); muForcingX3.DefineVar("x3",&muX3);
-
- muDeltaT = deltaT;
-
- muForcingX1.DefineVar("dt",&muDeltaT);
- muForcingX2.DefineVar("dt",&muDeltaT);
- muForcingX3.DefineVar("dt",&muDeltaT);
-
- muNu = (1.0/3.0)*(1.0/collFactor - 1.0/2.0);
-
- muForcingX1.DefineVar("nu",&muNu);
- muForcingX2.DefineVar("nu",&muNu);
- muForcingX3.DefineVar("nu",&muNu);
-
- LBMReal forcingX1 = 0;
- LBMReal forcingX2 = 0;
- LBMReal forcingX3 = 0;
- }*/
forcingX1 = 0.0;
forcingX2 = 0.0;
forcingX3 = 0.0;
/////////////////////////////////////
- localDistributionsF =
- dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
- nonLocalDistributionsF =
- dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
- zeroDistributionsF =
- dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
-
- localDistributionsH =
- dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getLocalDistributions();
- nonLocalDistributionsH =
- dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getNonLocalDistributions();
- zeroDistributionsH =
- dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getZeroDistributions();
+ localDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
+ nonLocalDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
+ zeroDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
- // phaseField = dataSet->getPhaseField();
+ localDistributionsH = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getLocalDistributions();
+ nonLocalDistributionsH = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getNonLocalDistributions();
+ zeroDistributionsH = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getZeroDistributions();
SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
@@ -142,19 +111,11 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
int maxX2 = bcArrayMaxX2 - ghostLayerWidth;
int maxX3 = bcArrayMaxX3 - ghostLayerWidth;
- //#pragma omp parallel num_threads(8)
- {
- // int i = omp_get_thread_num();
- // printf_s("Hello from thread %d\n", i);
- //}
- //#pragma omp for
-
CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField(
new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, -999.0));
CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr divU(
new CbArray3D<LBMReal, IndexerX3X2X1>(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3, 0.0));
- // CbArray3D<LBMReal> phaseField(bcArrayMaxX1, bcArrayMaxX2, bcArrayMaxX3,-999);
for (int x3 = 0; x3 <= maxX3; x3++) {
for (int x2 = 0; x2 <= maxX2; x2++) {
@@ -192,10 +153,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
LBMReal mfcca = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p);
LBMReal mfbbb = (*this->zeroDistributionsH)(x1, x2, x3);
- // LBMReal phase = h[REST] + h[E] + h[W] + h[N] + h[S] + h[T] + h[B] + h[NE] + h[SW] + h[SE] +
- // h[NW] + h[TE] + h[BW] + h[BE] + h[TW] + h[TN] + h[BS] + h[BN] + h[TS] + h[TNE] + h[TNW] +
- //h[TSE] + h[TSW] + h[BNE] + h[BNW] + h[BSE] + h[BSW]; if (phase > 1.0) phase = 1.0e0;
- //(*phaseField)(x1,x2,x3) = phase;
(*phaseField)(x1, x2, x3) = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca) +
(mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) +
(mfbaa + mfbac + mfbca + mfbcc) + (mfabb + mfcbb) +
@@ -207,31 +164,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
LBMReal collFactorM;
LBMReal forcingTerm[D3Q27System::ENDF + 1];
- // LBMReal m000, m100, m010, m001, m110, m101, m011, m200, m020, m002, m120, m102, m210, m012, m201, m021, m111,
- // m220, m202, m022, m211, m121, m112, m221, m212, m122, m222; LBMReal k000, k100, k010, k001, k110, k101, k011,
- // k200, k020, k002, k120, k102, k210, k012, k201, k021, k111, k220, k202, k022, k211, k121, k112, k221, k212,
- // k122, k222; LBMReal c000, c100, c010, c001, c110, c101, c011, c200, c020, c002, c120, c102, c210, c012, c201,
- // c021, c111, c220, c202, c022, c211, c121, c112, c221, c212, c122, c222;
-
- // LBMReal k200_pl_k020_pl_k002, k200_mi_k020, k200_mi_k002, k210_pl_k012, k210_mi_k012, k201_pl_k021,
- // k201_mi_k021, k120_pl_k102, k120_mi_k102, k220_pl_k202_pl_k022,
- // k220_mi2_k202_pl_k022, k220_pl_k202_mi2_k022;
-
- // LBMReal c200_pl_c020_pl_c002, c200_mi_c020, c200_mi_c002, c210_pl_c012, c210_mi_c012, c201_pl_c021,
- // c201_mi_c021, c120_pl_c102, c120_mi_c102, c220_pl_c202_pl_c022,
- // c220_mi2_c202_pl_c022, c220_pl_c202_mi2_c022;
-
- LBMReal w1, w2, w3, w4, w5, w6, w7, w8, w9, w10;
-
- w2 = 1.0;
- w3 = 1.0;
- w4 = 1.0;
- w5 = 1.0;
- w6 = 1.0;
- w7 = 1.0;
- w8 = 1.0;
- w9 = 1.0;
- w10 = 1.0;
for (int x3 = minX3; x3 < maxX3; x3++) {
for (int x2 = minX2; x2 < maxX2; x2++) {
@@ -262,35 +194,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
// a b c
//-1 0 1
- /*
- phi[REST] = (phaseField)(x1,x2,x3);
- phi[E ] = (phaseField)(x1 + DX1[E ], x2 + DX2[E ], x3 + DX3[E ]);
- phi[N ] = (phaseField)(x1 + DX1[N ], x2 + DX2[N ], x3 + DX3[N ]);
- phi[T ] = (phaseField)(x1 + DX1[T ], x2 + DX2[T ], x3 + DX3[T ]);
- phi[W ] = (phaseField)(x1 + DX1[W ], x2 + DX2[W ], x3 + DX3[W ]);
- phi[S ] = (phaseField)(x1 + DX1[S ], x2 + DX2[S ], x3 + DX3[S ]);
- phi[B ] = (phaseField)(x1 + DX1[B ], x2 + DX2[B ], x3 + DX3[B ]);
- phi[NE ] = (phaseField)(x1 + DX1[NE ], x2 + DX2[NE ], x3 + DX3[NE ]);
- phi[NW ] = (phaseField)(x1 + DX1[NW ], x2 + DX2[NW ], x3 + DX3[NW ]);
- phi[TE ] = (phaseField)(x1 + DX1[TE ], x2 + DX2[TE ], x3 + DX3[TE ]);
- phi[TW ] = (phaseField)(x1 + DX1[TW ], x2 + DX2[TW ], x3 + DX3[TW ]);
- phi[TN ] = (phaseField)(x1 + DX1[TN ], x2 + DX2[TN ], x3 + DX3[TN ]);
- phi[TS ] = (phaseField)(x1 + DX1[TS ], x2 + DX2[TS ], x3 + DX3[TS ]);
- phi[SW ] = (phaseField)(x1 + DX1[SW ], x2 + DX2[SW ], x3 + DX3[SW ]);
- phi[SE ] = (phaseField)(x1 + DX1[SE ], x2 + DX2[SE ], x3 + DX3[SE ]);
- phi[BW ] = (phaseField)(x1 + DX1[BW ], x2 + DX2[BW ], x3 + DX3[BW ]);
- phi[BE ] = (phaseField)(x1 + DX1[BE ], x2 + DX2[BE ], x3 + DX3[BE ]);
- phi[BS ] = (phaseField)(x1 + DX1[BS ], x2 + DX2[BS ], x3 + DX3[BS ]);
- phi[BN ] = (phaseField)(x1 + DX1[BN ], x2 + DX2[BN ], x3 + DX3[BN ]);
- phi[BSW] = (phaseField)(x1 + DX1[BSW], x2 + DX2[BSW], x3 + DX3[BSW]);
- phi[BSE] = (phaseField)(x1 + DX1[BSE], x2 + DX2[BSE], x3 + DX3[BSE]);
- phi[BNW] = (phaseField)(x1 + DX1[BNW], x2 + DX2[BNW], x3 + DX3[BNW]);
- phi[BNE] = (phaseField)(x1 + DX1[BNE], x2 + DX2[BNE], x3 + DX3[BNE]);
- phi[TNE] = (phaseField)(x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE]);
- phi[TNW] = (phaseField)(x1 + DX1[TNW], x2 + DX2[TNW], x3 + DX3[TNW]);
- phi[TSE] = (phaseField)(x1 + DX1[TSE], x2 + DX2[TSE], x3 + DX3[TSE]);
- phi[TSW] = (phaseField)(x1 + DX1[TSW], x2 + DX2[TSW], x3 + DX3[TSW]);
- */
findNeighbors(phaseField, x1, x2, x3);
LBMReal mfcbb = (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3);
@@ -325,62 +228,20 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
LBMReal rhoH = 1.0;
LBMReal rhoL = 1.0 / densityRatio;
- // LBMReal rhoToPhi = (1.0 - 1.0/densityRatio);
LBMReal rhoToPhi = (rhoH - rhoL) / (phiH - phiL);
- // collFactorM = phi[REST]*collFactorL + (1-phi[REST])*collFactorG;
- // collFactorM = phi[REST]*collFactorG + (1-phi[REST])*collFactorL;
-
- // LBMReal tauH = 1.0;
- // LBMReal di = sqrt(8*kappa/beta);
-
LBMReal dX1_phi = gradX1_phi();
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;
-
collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
- /*if ((phi[REST] > 0.1)||(phi[REST] < 0.9))
- {
- collFactorM*=(1.0-denom);
- }*/
-
- w1 = collFactorM;
-
- /*dX1_phi = -normX1*((phi[REST]>phiH || phi[REST]<phiL) ? 0.0 : 4*(phi[REST] - phiL)*(phi[REST]
- - phiH)/di); dX2_phi = -normX2*((phi[REST]>phiH || phi[REST]<phiL) ? 0.0 : 4*(phi[REST] -
- phiL)*(phi[REST] - phiH)/di); dX3_phi = -normX3*((phi[REST]>phiH || phi[REST]<phiL) ? 0.0 :
- 4*(phi[REST] - phiL)*(phi[REST] - phiH)/di);*/
- // UbTupleDouble3 coords = grid->getNodeCoordinates(block, x1, x2, x3);
- /*Block3D bl = this->block();
-
- int wX1 = bl->getX1() + x1;
- int wX2 = bl->getX2() + x2;
- int wX3 = bl->getX3() + x3;*/
-
- /*if (wX3 >= 30.0)
- {
- dX1_phi = 0.0;
- dX2_phi = 0.0;
- dX3_phi = 0.0;
- }*/
-
- LBMReal mu =
- 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
-
- // LBMReal rhoToPhi = (1.0/densityRatio - 1.0);
+ LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
//----------- Calculating Macroscopic Values -------------
-
- // LBMReal rho = phi[REST] + (1.0 - phi[REST])*1.0/densityRatio;
LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
- // LBMReal rho = phi[REST]*1.0/densityRatio + (1.0 - phi[REST]);
if (withForcing) {
// muX1 = static_cast<double>(x1-1+ix1*maxX1);
@@ -419,62 +280,19 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
(rho * c1o3) +
(mu * dX3_phi + forcingX3) / (2 * rho);
- // LBMReal p1 = (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 +
- //(ux*rhoToPhi*dX1_phi*c1o3 + uy*rhoToPhi*dX2_phi*c1o3 + uz*rhoToPhi*dX3_phi*c1o3)/2.0;
-
- // vvx = 0.0; vvy = 0.0; vvz = 0.0;
//--------------------------------------------------------
LBMReal ux2 = ux * ux;
LBMReal uy2 = uy * uy;
LBMReal uz2 = uz * uz;
- // LBMReal ux_uy = ux*uy;
- // LBMReal ux_uz = ux*uz;
- // LBMReal uy_uz = uy*uz;
- // LBMReal ux_uy_uz = ux*uy*uz;
-
- /*
- //----------- Calculating Forcing Terms -------------
- LBMReal forcingTerm1 = (ux*mu*dX1_phi + uy*mu*dX2_phi + uz*mu*dX3_phi);
- for (int dir = STARTF; dir < (ENDF+1); dir++)
- {
- if (dir != REST)
- {
- LBMReal dirGrad_phi = (phi[dir] - phi[INVDIR[dir]])/2.0;
- forcingTerm[dir] = (c1o3*rhoToPhi*dirGrad_phi +
- mu*dirGrad_phi)*(DX1[dir]*ux + DX2[dir]*uy + DX3[dir]*uz)*WEIGTH[dir]/c1o3 +
- mu*dirGrad_phi*WEIGTH[dir] - (forcingTerm1)*WEIGTH[dir];
- }
- else
- {
- forcingTerm[REST] = -(forcingTerm1)*WEIGTH[REST];
- }
- }
- //--------------------------------------------------------
- */
//----------- Calculating Forcing Terms * -------------
- // LBMReal forcingTerm1 = (ux*mu*dX1_phi + uy*mu*dX2_phi + uz*mu*dX3_phi);
for (int dir = STARTF; dir <= (FENDDIR); dir++) {
LBMReal velProd = DX1[dir] * ux + DX2[dir] * uy + DX3[dir] * uz;
LBMReal velSq1 = velProd * velProd;
LBMReal gamma = WEIGTH[dir] * (1.0 + 3 * velProd + 4.5 * velSq1 - 1.5 * (ux2 + uy2 + uz2));
- // forcingTerm[dir] = (DX1[dir] - ux)*((gamma - WEIGTH[dir])*c1o3*rhoToPhi*dX1_phi +
- // gamma*mu*dX1_phi) +
- // (DX2[dir] - uy)*((gamma - WEIGTH[dir])*c1o3*rhoToPhi*dX2_phi + gamma*mu*dX2_phi) +
- // (DX3[dir] - uz)*((gamma - WEIGTH[dir])*c1o3*rhoToPhi*dX3_phi + gamma*mu*dX3_phi);
-
LBMReal fac1 = (gamma - WEIGTH[dir]) * c1o3 * rhoToPhi;
- // LBMReal dirGrad_phi = (phi[dir] - phi[INVDIR[dir]])/2.0;
- // LBMReal dirGrad_phi = DX1[dir]*dX1_phi + DX2[dir]*dX2_phi + DX3[dir]*dX3_phi;
-
- /*forcingTerm[dir] = (- (ux)*(fac1*dX1_phi + gamma*mu*dX1_phi) -
- (uy)*(fac1*dX2_phi + gamma*mu*dX2_phi) -
- (uz)*(fac1*dX3_phi + gamma*mu*dX3_phi)) + (fac1*dirGrad_phi + gamma*mu*dirGrad_phi +
- DX1[dir]*forcingX1 + DX2[dir]*forcingX2 + DX3[dir]*forcingX3);*/
forcingTerm[dir] = ((-ux) * (fac1 * dX1_phi + gamma * (mu * dX1_phi + forcingX1)) +
(-uy) * (fac1 * dX2_phi + gamma * (mu * dX2_phi + forcingX2)) +
@@ -485,9 +303,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
}
LBMReal gamma = WEIGTH[REST] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
- /*forcingTerm[REST] = -(ux)*((gamma - WEIGTH[REST])*c1o3*rhoToPhi*dX1_phi + gamma*mu*dX1_phi) -
- (uy)*((gamma - WEIGTH[REST])*c1o3*rhoToPhi*dX2_phi + gamma*mu*dX2_phi) -
- (uz)*((gamma - WEIGTH[REST])*c1o3*rhoToPhi*dX3_phi + gamma*mu*dX3_phi);*/
LBMReal fac1 = (gamma - WEIGTH[REST]) * c1o3 * rhoToPhi;
forcingTerm[REST] = (-ux) * (fac1 * dX1_phi + gamma * (mu * dX1_phi + forcingX1)) +
(-uy) * (fac1 * dX2_phi + gamma * (mu * dX2_phi + forcingX2)) +
@@ -495,37 +310,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
//--------------------------------------------------------
- /*
- f1[E ] = (g[E ] + 0.5*forcingTerm[E ])/c1o3 - (p1 - rho*c1o3)*WEIGTH[E ]/c1o3;
- f1[N ] = (g[N ] + 0.5*forcingTerm[N ])/c1o3 - (p1 - rho*c1o3)*WEIGTH[N ]/c1o3;
- f1[T ] = (g[T ] + 0.5*forcingTerm[T ])/c1o3 - (p1 - rho*c1o3)*WEIGTH[T ]/c1o3;
- f1[NE ] = (g[NE ] + 0.5*forcingTerm[NE ])/c1o3 - (p1 - rho*c1o3)*WEIGTH[NE
- ]/c1o3; f1[NW ] = (g[NW ] + 0.5*forcingTerm[NW ])/c1o3 - (p1 - rho*c1o3)*WEIGTH[NW ]/c1o3;
- f1[TE ] = (g[TE ] + 0.5*forcingTerm[TE ])/c1o3 - (p1 - rho*c1o3)*WEIGTH[TE
- ]/c1o3; f1[TW ] = (g[TW ] + 0.5*forcingTerm[TW ])/c1o3 - (p1 - rho*c1o3)*WEIGTH[TW ]/c1o3;
- f1[TN ] = (g[TN ] + 0.5*forcingTerm[TN ])/c1o3 - (p1 - rho*c1o3)*WEIGTH[TN
- ]/c1o3; f1[TS ] = (g[TS ] + 0.5*forcingTerm[TS ])/c1o3 - (p1 - rho*c1o3)*WEIGTH[TS ]/c1o3;
- f1[TNE] = (g[TNE] + 0.5*forcingTerm[TNE])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[TNE]/c1o3; f1[TNW] = (g[TNW] + 0.5*forcingTerm[TNW])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[TNW]/c1o3; f1[TSE] = (g[TSE] + 0.5*forcingTerm[TSE])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[TSE]/c1o3; f1[TSW] = (g[TSW] + 0.5*forcingTerm[TSW])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[TSW]/c1o3; f1[W ] = (g[W ] + 0.5*forcingTerm[W ])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[W ]/c1o3; f1[S ] = (g[S ] + 0.5*forcingTerm[S ])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[S ]/c1o3; f1[B ] = (g[B ] + 0.5*forcingTerm[B ])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[B ]/c1o3; f1[SW ] = (g[SW ] + 0.5*forcingTerm[SW ])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[SW ]/c1o3; f1[SE ] = (g[SE ] + 0.5*forcingTerm[SE ])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[SE ]/c1o3; f1[BW ] = (g[BW ] + 0.5*forcingTerm[BW ])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[BW ]/c1o3; f1[BE ] = (g[BE ] + 0.5*forcingTerm[BE ])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[BE ]/c1o3; f1[BS ] = (g[BS ] + 0.5*forcingTerm[BS ])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[BS ]/c1o3; f1[BN ] = (g[BN ] + 0.5*forcingTerm[BN ])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[BN ]/c1o3; f1[BSW] = (g[BSW] + 0.5*forcingTerm[BSW])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[BSW]/c1o3; f1[BSE] = (g[BSE] + 0.5*forcingTerm[BSE])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[BSE]/c1o3; f1[BNW] = (g[BNW] + 0.5*forcingTerm[BNW])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[BNW]/c1o3; f1[BNE] = (g[BNE] + 0.5*forcingTerm[BNE])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[BNE]/c1o3; f1[REST] = (g[REST] + 0.5*forcingTerm[REST])/c1o3 - (p1 -
- rho*c1o3)*WEIGTH[REST]/c1o3;
- */
-
mfcbb = 3.0 * (mfcbb + 0.5 * forcingTerm[E]) / rho; //-(3.0*p1 - rho)*WEIGTH[E ];
mfbcb = 3.0 * (mfbcb + 0.5 * forcingTerm[N]) / rho; //-(3.0*p1 - rho)*WEIGTH[N ];
mfbbc = 3.0 * (mfbbc + 0.5 * forcingTerm[T]) / rho; //-(3.0*p1 - rho)*WEIGTH[T ];
@@ -559,25 +343,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
(mfbaa + mfbac + mfbca + mfbcc) + (mfabb + mfcbb) + (mfbab + mfbcb) +
(mfbba + mfbbc) + mfbbb;
- /*
- //forcing
- ///////////////////////////////////////////////////////////////////////////////////////////
- if (withForcing)
- {
- muX1 = static_cast<double>(x1-1+ix1*maxX1);
- muX2 = static_cast<double>(x2-1+ix2*maxX2);
- muX3 = static_cast<double>(x3-1+ix3*maxX3);
-
- forcingX1 = muForcingX1.Eval();
- forcingX2 = muForcingX2.Eval();
- forcingX3 = muForcingX3.Eval();
-
- vvx += forcingX1*deltaT*0.5; // X
- vvy += forcingX2*deltaT*0.5; // Y
- vvz += forcingX3*deltaT*0.5; // Z
- }
- ///////////////////////////////////////////////////////////////////////////////////////////
- */
LBMReal oMdrho, m0, m1, m2;
@@ -849,11 +614,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
LBMReal O6 = 1.;
// 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 LBMReal CUMbbc = mfbbc - ((mfaac + c1o3 * oMdrho) * mfbba + 2. *
- // mfbab * mfabb); // till 18.05.2015
-
LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + 2. * mfbba * mfbab);
LBMReal CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + 2. * mfbba * mfabb);
LBMReal CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + 2. * mfbab * mfabb);
@@ -865,12 +625,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
LBMReal CUMacc = mfacc - ((mfaac * mfaca + 2. * mfabb * mfabb) +
c1o3 * (mfaac + mfaca) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
- // LBMReal CUMcca = mfcca - ((mfcaa * mfaca + 2. * mfbba * mfbba) + c1o3 * (mfcaa + mfaca) *
- // oMdrho + c1o9*(-p1/c1o3)*oMdrho); LBMReal CUMcac = mfcac - ((mfcaa * mfaac + 2. * mfbab *
- // mfbab) + c1o3 * (mfcaa + mfaac) * oMdrho + c1o9*(-p1/c1o3)*oMdrho); LBMReal CUMacc = mfacc -
- // ((mfaac * mfaca + 2. * mfabb * mfabb) + c1o3 * (mfaac + mfaca) * oMdrho +
- // c1o9*(-p1/c1o3)*oMdrho);
-
// Cum 5.
LBMReal CUMbcc = mfbcc -
(mfaac * mfbca + mfaca * mfbac + 4. * mfabb * mfbbb +
@@ -911,14 +665,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
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;
-
- LBMReal strainMag = sqrt(2*(dxux*dxux + dyuy*dyuy + dzuz*dzuz) + Dxy*Dxy + Dxz*Dxz + Dyz*Dyz);
- LBMReal intVis = 3*abs(denom - 1e-9)*strainMag;
- LBMReal fluidVis = (1.0/collFactorM - 0.5)/3.0;
- collFactorM = 1.0/((fluidVis + intVis)*3.0 + 0.5);*/
(*divU)(x1, x2, x3) = dxux + dyuy + dzuz;
// relax
@@ -990,10 +736,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
// back cumulants to central moments
// 4.
- // mfcbb = CUMcbb + ((mfcaa + c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
- // mfbcb = CUMbcb + ((mfaca + c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015
- // mfbbc = CUMbbc + ((mfaac + c1o3 * oMdrho) * mfbba + 2. * mfbab * mfabb); // till 18.05.2015
-
mfcbb = CUMcbb + ((mfcaa + c1o3) * mfabb + 2. * mfbba * mfbab);
mfbcb = CUMbcb + ((mfaca + c1o3) * mfbab + 2. * mfbba * mfabb);
mfbbc = CUMbbc + ((mfaac + c1o3) * mfbba + 2. * mfbab * mfabb);
@@ -1005,11 +747,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
mfacc = CUMacc + (mfaac * mfaca + 2. * mfabb * mfabb) + c1o3 * (mfaac + mfaca) * oMdrho +
c1o9 * (oMdrho - 1) * oMdrho;
- // mfcca = CUMcca + (mfcaa * mfaca + 2. * mfbba * mfbba) + c1o3 * (mfcaa + mfaca) * oMdrho +
- // c1o9*(-p1/c1o3)*oMdrho; mfcac = CUMcac + (mfcaa * mfaac + 2. * mfbab * mfbab) + c1o3 * (mfcaa
- // + mfaac) * oMdrho + c1o9*(-p1/c1o3)*oMdrho; mfacc = CUMacc + (mfaac * mfaca + 2. * mfabb *
- // mfabb) + c1o3 * (mfaac + mfaca) * oMdrho + c1o9*(-p1/c1o3)*oMdrho;
-
// 5.
mfbcc = CUMbcc +
(mfaac * mfbca + mfaca * mfbac + 4. * mfabb * mfbbb +
@@ -1266,10 +1003,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
(mfbaa + mfbac + mfbca + mfbcc) + (mfabb + mfcbb) + (mfbab + mfbcb) +
(mfbba + mfbbc) + mfbbb;
- /*LBMReal rho_post = f1[REST] + f1[E] + f1[W] + f1[N] + f1[S] + f1[T] + f1[B] + f1[NE] + f1[SW]
- + f1[SE] + f1[NW] + f1[TE] + f1[BW] + f1[BE] + f1[TW] + f1[TN] + f1[BS] + f1[BN] + f1[TS] +
- f1[TNE] + f1[TNW] + f1[TSE] + f1[TSW] + f1[BNE] + f1[BNW] + f1[BSE] + f1[BSW]; */
- // LBMReal dif = fabs(rho - rho_post);
LBMReal dif = rho1 - rho_post;
#ifdef SINGLEPRECISION
if (dif > 10.0E-7 || dif < -10.0E-7)
@@ -1282,9 +1015,6 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
UbSystem::toString(rho_post) + " dif=" + UbSystem::toString(dif) +
" rho is not correct for node " + UbSystem::toString(x1) + "," +
UbSystem::toString(x2) + "," + UbSystem::toString(x3)));
- // UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): rho is not correct for node
- // "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3));
- // exit(EXIT_FAILURE);
}
#endif
@@ -1385,41 +1115,22 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
h[REST] = (*this->zeroDistributionsH)(x1, x2, x3);
- // LBMReal denom = sqrt(dX1_phi*dX1_phi + dX2_phi*dX2_phi + dX3_phi*dX3_phi) + 1e-15;
- // LBMReal di = sqrt(8*kappa/beta);
- // LBMReal tauH1 = 3.0*mob + 0.5;
for (int dir = STARTF; dir < (ENDF + 1); dir++) {
LBMReal velProd = DX1[dir] * ux + DX2[dir] * uy + DX3[dir] * uz;
LBMReal velSq1 = velProd * velProd;
LBMReal hEq; //, gEq;
if (dir != REST) {
- // LBMReal dirGrad_phi = DX1[dir]*dX1_phi+DX2[dir]*dX2_phi+DX3[dir]*dX3_phi;
LBMReal dirGrad_phi = (phi[dir] - phi[INVDIR[dir]]) / 2.0;
- LBMReal hSource = (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST]) * (dirGrad_phi) /
- denom; // + phi[REST]*(dxux + dyuy + dzuz);
-
- // LBMReal hSource =((phi[REST]>phiH || phi[REST]<phiL) ? 0.1 : 1.0)
- // * 3.0*mob*(-4.0)/di*(phi[REST] - phiL)*(phi[REST] - phiH)*(dirGrad_phi)/denom; LBMReal
- // hSource = 3.0*mob*(-4.0)/di*(phi[REST] - phiL)*(phi[REST] - phiH)*(dirGrad_phi)/denom;
- hEq = phi[REST] * WEIGTH[dir] *
- (1.0 + 3.0 * velProd + 4.5 * velSq1 - 1.5 * (ux2 + uy2 + uz2)) +
- hSource * WEIGTH[dir];
- // gEq = rho*WEIGTH[dir]*(1 + 3*velProd + 4.5*velSq1 - 1.5*(vx2+vy2+vz2))*c1o3 +
- // (p1-rho*c1o3)*WEIGTH[dir]; h[dir] = hEq; //h[dir] - (h[dir] - hEq)/(tauH + 0.5)); ///
+ LBMReal hSource = (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST]) * (dirGrad_phi) / denom;
+ hEq = phi[REST] * WEIGTH[dir] * (1.0 + 3.0 * velProd + 4.5 * velSq1 - 1.5 * (ux2 + uy2 + uz2)) + hSource * WEIGTH[dir];
+
// This corresponds with the collision factor of 1.0 which equals (tauH + 0.5).
- h[dir] =
- h[dir] - (h[dir] - hEq) / (tauH); // + WEIGTH[dir]*phi[REST]*(dxux + dyuy + dzuz);
- // h[dir] = h[dir] - (h[dir] - hEq)/(tauH1);
- // g[dir] = g[dir] - collFactorM*(g[dir]-gEq) + 0.5*forcingTerm[dir];
+ h[dir] = h[dir] - (h[dir] - hEq) / (tauH);
} else {
hEq = phi[REST] * WEIGTH[REST] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
- // gEq = rho*WEIGTH[dir]*(1 + 3*velProd + 4.5*velSq1 - 1.5*(vx2+vy2+vz2))*c1o3 +
- // (p1-rho*c1o3)*WEIGTH[dir]; h[dir] = hEq;
- h[REST] = h[REST] -
- (h[REST] - hEq) / (tauH); // + WEIGTH[REST]*phi[REST]*(dxux + dyuy + dzuz);
- // g[dir] = g[dir] - collFactorM*(g[dir]-gEq) + 0.5*forcingTerm[dir];
+ h[REST] = h[REST] - (h[REST] - hEq) / (tauH);
}
}
@@ -1460,7 +1171,7 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
}
dataSet->setPhaseField(divU);
}
-}
+
//////////////////////////////////////////////////////////////////////////
LBMReal MultiphaseCumulantLBMKernel::gradX1_phi()
@@ -1502,17 +1213,12 @@ LBMReal MultiphaseCumulantLBMKernel::nabla2_phi()
}
return 6.0 * sum;
}
-///// Commnets neeeded ////////
void MultiphaseCumulantLBMKernel::computePhasefield()
{
using namespace D3Q27System;
SPtr<DistributionArray3D> distributionsH = dataSet->getHdistributions();
- // const int bcArrayMaxX1 = (int)distributionsH->getNX1();
- // const int bcArrayMaxX2 = (int)distributionsH->getNX2();
- // const int bcArrayMaxX3 = (int)distributionsH->getNX3();
-
int minX1 = ghostLayerWidth;
int minX2 = ghostLayerWidth;
int minX3 = ghostLayerWidth;
@@ -1559,61 +1265,10 @@ void MultiphaseCumulantLBMKernel::computePhasefield()
h[BNE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p);
h[REST] = (*this->zeroDistributionsH)(x1, x2, x3);
-
- /*(*this->phaseField)(x1,x2,x3) = h[REST] + h[E] + h[W] + h[N] + h[S] + h[T] + h[B] + h[NE] + h[SW]
- + h[SE] + h[NW] + h[TE] + h[BW] +
- h[BE] + h[TW] + h[TN] + h[BS] + h[BN] + h[TS] + h[TNE] + h[TNW] + h[TSE] + h[TSW] + h[BNE] +
- h[BNW] + h[BSE] + h[BSW];*/
}
}
}
}
- //----------------------------------------------------------
-
- /*
- /////// Filling ghost nodes for FD computations //////////
- for(int x1 = minX1; x1 < maxX1; x1++)
- {
- for(int x2 = minX2; x2 < maxX2; x2++)
- {
- int x3 = 0;
- (*phaseField)(x1, x2, x3) = (*phaseField)(x1, x2, maxX3-1);
- x3 = maxX3;
- (*phaseField)(x1, x2, x3) = (*phaseField)(x1, x2, minX3);
- }
- }
- for(int x2 = minX2; x2 < maxX2; x2++)
- {
- for(int x3 = minX3; x3 < maxX3; x3++)
- {
- int x1 = 0;
- (*phaseField)(x1, x2, x3) = (*phaseField)(maxX1-1, x2, x3);
- x1 = maxX1;
- (*phaseField)(x1, x2, x3) = (*phaseField)(minX1, x2, x3);
- }
- }
- for(int x1 = minX1; x1 < maxX1; x1++)
- {
- for(int x3 = minX3; x3 < maxX3; x3++)
- {
- int x2 = 0;
- (*phaseField)(x1, x2, x3) = (*phaseField)(x1, maxX2-1, x3);
- x2 = maxX2;
- (*phaseField)(x1, x2, x3) = (*phaseField)(x1, minX2, x3);
- }
- }
- (*phaseField)(0, 0, 0 ) = (*phaseField)(maxX1-1, maxX2-1, maxX3-1);
- (*phaseField)(0, 0, maxX3) = (*phaseField)(maxX1-1, maxX2-1, minX3 );
- (*phaseField)(0, maxX2, 0 ) = (*phaseField)(maxX1-1, minX2, maxX3-1 );
- (*phaseField)(0, maxX2, maxX3) = (*phaseField)(maxX1-1, minX2, minX3 );
-
- (*phaseField)(maxX1, 0, 0 ) = (*phaseField)(minX1, maxX2-1, maxX3-1);
- (*phaseField)(maxX1, 0, maxX3) = (*phaseField)(minX1, maxX2-1, minX3 );
- (*phaseField)(maxX1, maxX2, 0 ) = (*phaseField)(minX1, minX2, maxX3-1 );
- (*phaseField)(maxX1, maxX2, maxX3) = (*phaseField)(minX1, minX2, minX3 );
-
- /////////////////////////////////////////////////////////
- */
}
void MultiphaseCumulantLBMKernel::findNeighbors(CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr ph, int x1, int x2,
@@ -1625,124 +1280,17 @@ void MultiphaseCumulantLBMKernel::findNeighbors(CbArray3D<LBMReal, IndexerX3X2X1
phi[REST] = (*ph)(x1, x2, x3);
- // LBMReal a = -0.5*sqrt(2*beta/kappa)*cos(contactAngle*UbMath::PI/180);
- // LBMReal a1 = 1 + a;
-
for (int k = FSTARTDIR; k <= FENDDIR; k++) {
if (!bcArray->isSolid(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k])) {
phi[k] = (*ph)(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k]);
} else {
- /*
- if (phi[REST] < 1e-2)
- {
- phi[k] = (*ph)(x1 + DX1[INVDIR[k]], x2 + DX2[INVDIR[k]], x3 + DX3[INVDIR[k]]);
- }
- else
- {
- LBMReal phi_f = (*ph)(x1 + DX1[k], x2, x3 + DX3[k]);
- phi[k] = (a1 - sqrt(a1*a1 - 4*a*phi_f) )/a - phi_f;
- }
- */
-
- phi[k] = (*ph)(x1, x2, x3);
-
- // if (bcArray->isSolid(x1 + DX1[k], x2, x3))
- //{
- // phi[k] = (*ph)(x1, x2, x3);
- // //if (!bcArray->isSolid(x1 , x2 + DX2[k], x3 + DX3[k]))
- // //{
- // // //phi[k] = (*ph)(x1 , x2 + DX2[k], x3 + DX3[k]);
- // // LBMReal phi_f = (*ph)(x1 , x2 + DX2[k], x3 + DX3[k]);
- // // phi[k] = (a1 - sqrt(a1*a1 - 4*a*phi_f) )/a - phi_f;
- // //}
- // //else
- // //{
- // // phi[k] = (*ph)(x1, x2, x3);
- // //}
- //}
- //
- // if (bcArray->isSolid(x1 , x2 + DX2[k], x3))
- //{
- // phi[k] = (*ph)(x1, x2, x3);
- // //if (!bcArray->isSolid(x1 + DX1[k], x2 , x3 + DX3[k]))
- // //{
- // // //phi[k] = (*ph)(x1 + DX1[k], x2 , x3 + DX3[k]);
- // // LBMReal phi_f = (*ph)(x1 + DX1[k], x2 , x3 + DX3[k]);
- // // phi[k] = (a1 - sqrt(a1*a1 - 4*a*phi_f) )/a - phi_f;
- // //}
- // //else
- // //{
- // // phi[k] = (*ph)(x1, x2, x3);
- // //}
- //}
-
- // if (bcArray->isSolid(x1 , x2, x3+ DX3[k]))
- //{
- // if (!bcArray->isSolid(x1 + DX1[k], x2 + DX2[k], x3))
- // {
- // //phi[k] = (*ph)(x1 + DX1[k], x2 + DX2[k], x3);
- // LBMReal phi_f = (*ph)(x1 + DX1[k], x2 + DX2[k], x3);
- // phi[k] = (a1 - sqrt(a1*a1 - 4*a*phi_f) )/a - phi_f;
- // }
- // else
- // {
- // phi[k] = (*ph)(x1, x2, x3);
- // }
- //}
-
- /*if (bcArray->isSolid(x1 + DX1[k], x2, x3)) phi[k] = (*ph)(x1 , x2 + DX2[k], x3 + DX3[k]);
- if (bcArray->isSolid(x1 , x2 + DX2[k], x3)) phi[k] = (*ph)(x1 + DX1[k], x2 , x3 + DX3[k]);
- if (bcArray->isSolid(x1 , x2, x3+ DX3[k])) phi[k] = (*ph)(x1 + DX1[k], x2 + DX2[k], x3 );*/
-
- /*if (phi[REST] < 0.00001)
- {
- phi[k] = 0.0;
- }
- else
- {
- phi[k] = 0.5;
- }*/
-
- // phi[k] = 0.5;
- // phi[k] = (*ph)(x1, x2, x3);
- // phi[k] = (*ph)(x1 + DX1[INVDIR[k]], x2 + DX2[INVDIR[k]], x3 + DX3[INVDIR[k]]);
- }
+ }
}
-
- /*
- phi[E ] = (*ph)(x1 + DX1[E ], x2 + DX2[E ], x3 + DX3[E ]);
- phi[N ] = (*ph)(x1 + DX1[N ], x2 + DX2[N ], x3 + DX3[N ]);
- phi[T ] = (*ph)(x1 + DX1[T ], x2 + DX2[T ], x3 + DX3[T ]);
- phi[W ] = (*ph)(x1 + DX1[W ], x2 + DX2[W ], x3 + DX3[W ]);
- phi[S ] = (*ph)(x1 + DX1[S ], x2 + DX2[S ], x3 + DX3[S ]);
- phi[B ] = (*ph)(x1 + DX1[B ], x2 + DX2[B ], x3 + DX3[B ]);
- phi[NE ] = (*ph)(x1 + DX1[NE ], x2 + DX2[NE ], x3 + DX3[NE ]);
- phi[NW ] = (*ph)(x1 + DX1[NW ], x2 + DX2[NW ], x3 + DX3[NW ]);
- phi[TE ] = (*ph)(x1 + DX1[TE ], x2 + DX2[TE ], x3 + DX3[TE ]);
- phi[TW ] = (*ph)(x1 + DX1[TW ], x2 + DX2[TW ], x3 + DX3[TW ]);
- phi[TN ] = (*ph)(x1 + DX1[TN ], x2 + DX2[TN ], x3 + DX3[TN ]);
- phi[TS ] = (*ph)(x1 + DX1[TS ], x2 + DX2[TS ], x3 + DX3[TS ]);
- phi[SW ] = (*ph)(x1 + DX1[SW ], x2 + DX2[SW ], x3 + DX3[SW ]);
- phi[SE ] = (*ph)(x1 + DX1[SE ], x2 + DX2[SE ], x3 + DX3[SE ]);
- phi[BW ] = (*ph)(x1 + DX1[BW ], x2 + DX2[BW ], x3 + DX3[BW ]);
- phi[BE ] = (*ph)(x1 + DX1[BE ], x2 + DX2[BE ], x3 + DX3[BE ]);
- phi[BS ] = (*ph)(x1 + DX1[BS ], x2 + DX2[BS ], x3 + DX3[BS ]);
- phi[BN ] = (*ph)(x1 + DX1[BN ], x2 + DX2[BN ], x3 + DX3[BN ]);
- phi[BSW] = (*ph)(x1 + DX1[BSW], x2 + DX2[BSW], x3 + DX3[BSW]);
- phi[BSE] = (*ph)(x1 + DX1[BSE], x2 + DX2[BSE], x3 + DX3[BSE]);
- phi[BNW] = (*ph)(x1 + DX1[BNW], x2 + DX2[BNW], x3 + DX3[BNW]);
- phi[BNE] = (*ph)(x1 + DX1[BNE], x2 + DX2[BNE], x3 + DX3[BNE]);
- phi[TNE] = (*ph)(x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE]);
- phi[TNW] = (*ph)(x1 + DX1[TNW], x2 + DX2[TNW], x3 + DX3[TNW]);
- phi[TSE] = (*ph)(x1 + DX1[TSE], x2 + DX2[TSE], x3 + DX3[TSE]);
- phi[TSW] = (*ph)(x1 + DX1[TSW], x2 + DX2[TSW], x3 + DX3[TSW]);
- */
}
void MultiphaseCumulantLBMKernel::swapDistributions()
{
LBMKernel::swapDistributions();
dataSet->getHdistributions()->swap();
- // computePhasefield();
}
\ No newline at end of file