Skip to content
Snippets Groups Projects
Commit d62cb39e authored by kutscher's avatar kutscher
Browse files

clean multiphase kernal

parent ef020a02
No related branches found
No related tags found
1 merge request!25remove warnings in multiphase classes
Hide whitespace changes
Inline Side-by-side
src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
+ 14
− 466
View file @ d62cb39e
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment