From cdc71b302654b958b043ca2a0523bf2a38ce8547 Mon Sep 17 00:00:00 2001
From: "AMATERASU\\geier" <geier@irmb.tu-bs.de>
Date: Tue, 13 Apr 2021 14:17:41 +0200
Subject: [PATCH] Classical rescaling of phase field distribution to 1/3 rho as
common in Lee-literature. The source terms were removed from the moments and
returned to the distributions based on the original code of Safari.
---
.../MultiphaseVelocityBCAlgorithm.cpp | 8 +-
.../MultiphaseScratchCumulantLBMKernel.cpp | 439 +++++++++++-------
2 files changed, 277 insertions(+), 170 deletions(-)
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp
index 609f71308..58c359887 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp
@@ -104,15 +104,15 @@ void MultiphaseVelocityBCAlgorithm::applyBC()
phiBC = bcPtr->getBoundaryPhaseField();
- //D3Q27System::calcMultiphaseHeq(htemp, phiBC, vx1, vx2, vx3);
- D3Q27System::calcMultiphaseHeq(htemp, phiBC, bcPtr->getBoundaryVelocityX1(), bcPtr->getBoundaryVelocityX2(), bcPtr->getBoundaryVelocityX2());//30.03.2021 EQ phase field BC!
+ D3Q27System::calcMultiphaseHeq(htemp, phiBC, vx1, vx2, vx3);
+ //D3Q27System::calcMultiphaseHeq(htemp, phiBC, bcPtr->getBoundaryVelocityX1(), bcPtr->getBoundaryVelocityX2(), bcPtr->getBoundaryVelocityX2());//30.03.2021 EQ phase field BC!
for (int fdir = D3Q27System::STARTF; fdir<=D3Q27System::ENDF; fdir++)
{
if (bcPtr->hasVelocityBoundaryFlag(fdir))
{
- // LBMReal hReturn = htemp[fdir]+h[fdir]-heq[fdir];
+ LBMReal hReturn = htemp[fdir]+h[fdir]-heq[fdir];
//17.03.2021 Let us just set the plain eq
- LBMReal hReturn = htemp[fdir];
+ //LBMReal hReturn = htemp[fdir];
distributionsH->setDistributionForDirection(hReturn, nx1, nx2, nx3, fdir);
if (distributionsH2)
distributionsH2->setDistributionForDirection(hReturn, nx1, nx2, nx3, fdir);
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
index f0cfacbf8..ce8b66cb5 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
@@ -305,7 +305,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
///!test
collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
-
+ //collFactorM = phi[REST] - phiL < (phiH - phiL) * 0.05 ? collFactorG : collFactorL;
LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
@@ -328,35 +328,35 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
////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;
+ 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;
- 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;
+ 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;
+ mfbbb = (*this->zeroDistributionsF)(x1, x2, x3) / rho * c3;
@@ -384,6 +384,70 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
vvx += mu * dX1_phi*c1o2;
vvy += mu * dX2_phi * c1o2;
vvz += mu * dX3_phi * c1o2;
+
+
+ ///----Classic source term 8.4.2021
+
+ LBMReal ux2;
+ LBMReal uy2;
+ LBMReal uz2;
+ ux2 = vvx * vvx;
+ uy2 = vvy * vvy;
+ uz2 = vvz * vvz;
+ LBMReal forcingTerm[D3Q27System::ENDF + 1];
+ for (int dir = STARTF; dir <= (FENDDIR); dir++) {
+ LBMReal velProd = DX1[dir] * vvx + DX2[dir] * vvy + DX3[dir] * vvz;
+ 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] =
+ (-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 );
+ }
+
+ 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]
+
+ //--------------------------------------------------------
+
+
//forcing
///////////////////////////////////////////////////////////////////////////////////////////
@@ -687,12 +751,12 @@ 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 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 B = (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.
@@ -737,8 +801,8 @@ 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 = 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
+ 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;
@@ -783,7 +847,7 @@ 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 / 6.0) * (dX1_phi * vvx + dX2_phi * vvy + dX3_phi * vvz) * correctionScaling; // As in Hesam's code
+ 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;
@@ -848,9 +912,9 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
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 * B + (one - O4) * (CUMbbc);
- CUMbcb = -O4 * (one / collFactorM - c1o2) * Dxz * c1o3 * B + (one - O4) * (CUMbcb);
- CUMcbb = -O4 * (one / collFactorM - c1o2) * Dyz * c1o3 * B + (one - O4) * (CUMcbb);
+ 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);
@@ -1159,57 +1223,93 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
//////////////////////////////////////////////////////////////////////////
//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;
+
+ /////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]
+
+
+
+ ////////////////////
+
+
+ (*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
@@ -2349,88 +2449,95 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
- ///////////////////// PHASE-FIELD BGK SOLVER ///////////////////////////////
-
- //h[E] = (*this->localDistributionsH)(D3Q27System::ET_E, x1, x2, x3);
- //h[N] = (*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3);
- //h[T] = (*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3);
- //h[NE] = (*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3);
- //h[NW] = (*this->localDistributionsH)(D3Q27System::ET_NW, x1p, x2, x3);
- //h[TE] = (*this->localDistributionsH)(D3Q27System::ET_TE, x1, x2, x3);
- //h[TW] = (*this->localDistributionsH)(D3Q27System::ET_TW, x1p, x2, x3);
- //h[TN] = (*this->localDistributionsH)(D3Q27System::ET_TN, x1, x2, x3);
- //h[TS] = (*this->localDistributionsH)(D3Q27System::ET_TS, x1, x2p, x3);
- //h[TNE] = (*this->localDistributionsH)(D3Q27System::ET_TNE, x1, x2, x3);
- //h[TNW] = (*this->localDistributionsH)(D3Q27System::ET_TNW, x1p, x2, x3);
- //h[TSE] = (*this->localDistributionsH)(D3Q27System::ET_TSE, x1, x2p, x3);
- //h[TSW] = (*this->localDistributionsH)(D3Q27System::ET_TSW, x1p, x2p, x3);
-
- //h[W] = (*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3);
- //h[S] = (*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3);
- //h[B] = (*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p);
- //h[SW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_SW, x1p, x2p, x3);
- //h[SE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_SE, x1, x2p, x3);
- //h[BW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BW, x1p, x2, x3p);
- //h[BE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BE, x1, x2, x3p);
- //h[BS] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BS, x1, x2p, x3p);
- //h[BN] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BN, x1, x2, x3p);
- //h[BSW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BSW, x1p, x2p, x3p);
- //h[BSE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BSE, x1, x2p, x3p);
- //h[BNW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNW, x1p, x2, x3p);
- //h[BNE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p);
-
- //h[REST] = (*this->zeroDistributionsH)(x1, x2, x3);
-
- //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 = (phi[dir] - phi[INVDIR[dir]]) / 2.0;
- // 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);
-
- // } else {
- // hEq = phi[REST] * WEIGTH[REST] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
- // h[REST] = h[REST] - (h[REST] - hEq) / (tauH);
- // }
- //}
-
- //(*this->localDistributionsH)(D3Q27System::ET_E, x1, x2, x3) = h[D3Q27System::INV_E];
- //(*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3) = h[D3Q27System::INV_N];
- //(*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3) = h[D3Q27System::INV_T];
- //(*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3) = h[D3Q27System::INV_NE];
- //(*this->localDistributionsH)(D3Q27System::ET_NW, x1p, x2, x3) = h[D3Q27System::INV_NW];
- //(*this->localDistributionsH)(D3Q27System::ET_TE, x1, x2, x3) = h[D3Q27System::INV_TE];
- //(*this->localDistributionsH)(D3Q27System::ET_TW, x1p, x2, x3) = h[D3Q27System::INV_TW];
- //(*this->localDistributionsH)(D3Q27System::ET_TN, x1, x2, x3) = h[D3Q27System::INV_TN];
- //(*this->localDistributionsH)(D3Q27System::ET_TS, x1, x2p, x3) = h[D3Q27System::INV_TS];
- //(*this->localDistributionsH)(D3Q27System::ET_TNE, x1, x2, x3) = h[D3Q27System::INV_TNE];
- //(*this->localDistributionsH)(D3Q27System::ET_TNW, x1p, x2, x3) = h[D3Q27System::INV_TNW];
- //(*this->localDistributionsH)(D3Q27System::ET_TSE, x1, x2p, x3) = h[D3Q27System::INV_TSE];
- //(*this->localDistributionsH)(D3Q27System::ET_TSW, x1p, x2p, x3) = h[D3Q27System::INV_TSW];
-
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3) = h[D3Q27System::INV_W];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3) = h[D3Q27System::INV_S];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p) = h[D3Q27System::INV_B];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_SW, x1p, x2p, x3) = h[D3Q27System::INV_SW];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_SE, x1, x2p, x3) = h[D3Q27System::INV_SE];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_BW, x1p, x2, x3p) = h[D3Q27System::INV_BW];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_BE, x1, x2, x3p) = h[D3Q27System::INV_BE];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_BS, x1, x2p, x3p) = h[D3Q27System::INV_BS];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_BN, x1, x2, x3p) = h[D3Q27System::INV_BN];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_BSW, x1p, x2p, x3p) = h[D3Q27System::INV_BSW];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_BSE, x1, x2p, x3p) = h[D3Q27System::INV_BSE];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_BNW, x1p, x2, x3p) = h[D3Q27System::INV_BNW];
- //(*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p) = h[D3Q27System::INV_BNE];
-
- //(*this->zeroDistributionsH)(x1, x2, x3) = h[D3Q27System::REST];
-
- ///////////////////// END OF OLD BGK SOLVER ///////////////////////////////
+ /////////////////// PHASE-FIELD BGK SOLVER ///////////////////////////////
+//using namespace D3Q27System;
+
+ // h[E] = (*this->localDistributionsH)(D3Q27System::ET_E, x1, x2, x3);
+ // h[N] = (*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3);
+ // h[T] = (*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3);
+ // h[NE] = (*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3);
+ // h[NW] = (*this->localDistributionsH)(D3Q27System::ET_NW, x1p, x2, x3);
+ // h[TE] = (*this->localDistributionsH)(D3Q27System::ET_TE, x1, x2, x3);
+ // h[TW] = (*this->localDistributionsH)(D3Q27System::ET_TW, x1p, x2, x3);
+ // h[TN] = (*this->localDistributionsH)(D3Q27System::ET_TN, x1, x2, x3);
+ // h[TS] = (*this->localDistributionsH)(D3Q27System::ET_TS, x1, x2p, x3);
+ // h[TNE] = (*this->localDistributionsH)(D3Q27System::ET_TNE, x1, x2, x3);
+ // h[TNW] = (*this->localDistributionsH)(D3Q27System::ET_TNW, x1p, x2, x3);
+ // h[TSE] = (*this->localDistributionsH)(D3Q27System::ET_TSE, x1, x2p, x3);
+ // h[TSW] = (*this->localDistributionsH)(D3Q27System::ET_TSW, x1p, x2p, x3);
+
+ // h[W] = (*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3);
+ // h[S] = (*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3);
+ // h[B] = (*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p);
+ // h[SW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_SW, x1p, x2p, x3);
+ // h[SE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_SE, x1, x2p, x3);
+ // h[BW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BW, x1p, x2, x3p);
+ // h[BE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BE, x1, x2, x3p);
+ // h[BS] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BS, x1, x2p, x3p);
+ // h[BN] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BN, x1, x2, x3p);
+ // h[BSW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BSW, x1p, x2p, x3p);
+ // h[BSE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BSE, x1, x2p, x3p);
+ // h[BNW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNW, x1p, x2, x3p);
+ // h[BNE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p);
+
+ // h[REST] = (*this->zeroDistributionsH)(x1, x2, x3);
+ ////vvx *= 3;
+ ////vvy *= 3;
+ ////vvz *= 3;
+ ////vx2 = vvx * vvx;
+ ////vy2 = vvy * vvy;
+ ////vz2 = vvz * vvz;
+
+ // for (int dir = STARTF; dir < (ENDF + 1); dir++) {
+ // LBMReal velProd = DX1[dir] * vvx + DX2[dir] * vvy + DX3[dir] * vvz;
+ // LBMReal velSq1 = velProd * velProd;
+ // LBMReal hEq; //, gEq;
+
+ // if (dir != REST) {
+ // LBMReal dirGrad_phi = (phi[dir] - phi[INVDIR[dir]]) / 2.0;
+ // 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 * (vx2 + vy2 + vz2)) + 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);
+
+ // } else {
+ // hEq = phi[REST] * WEIGTH[REST] * (1.0 - 1.5 * (vx2 + vy2 + vz2));
+ // h[REST] = h[REST] - (h[REST] - hEq) / (tauH);
+ // }
+ // }
+
+ // (*this->localDistributionsH)(D3Q27System::ET_E, x1, x2, x3) = h[D3Q27System::INV_E];
+ // (*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3) = h[D3Q27System::INV_N];
+ // (*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3) = h[D3Q27System::INV_T];
+ // (*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3) = h[D3Q27System::INV_NE];
+ // (*this->localDistributionsH)(D3Q27System::ET_NW, x1p, x2, x3) = h[D3Q27System::INV_NW];
+ // (*this->localDistributionsH)(D3Q27System::ET_TE, x1, x2, x3) = h[D3Q27System::INV_TE];
+ // (*this->localDistributionsH)(D3Q27System::ET_TW, x1p, x2, x3) = h[D3Q27System::INV_TW];
+ // (*this->localDistributionsH)(D3Q27System::ET_TN, x1, x2, x3) = h[D3Q27System::INV_TN];
+ // (*this->localDistributionsH)(D3Q27System::ET_TS, x1, x2p, x3) = h[D3Q27System::INV_TS];
+ // (*this->localDistributionsH)(D3Q27System::ET_TNE, x1, x2, x3) = h[D3Q27System::INV_TNE];
+ // (*this->localDistributionsH)(D3Q27System::ET_TNW, x1p, x2, x3) = h[D3Q27System::INV_TNW];
+ // (*this->localDistributionsH)(D3Q27System::ET_TSE, x1, x2p, x3) = h[D3Q27System::INV_TSE];
+ // (*this->localDistributionsH)(D3Q27System::ET_TSW, x1p, x2p, x3) = h[D3Q27System::INV_TSW];
+
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3) = h[D3Q27System::INV_W];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3) = h[D3Q27System::INV_S];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p) = h[D3Q27System::INV_B];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_SW, x1p, x2p, x3) = h[D3Q27System::INV_SW];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_SE, x1, x2p, x3) = h[D3Q27System::INV_SE];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_BW, x1p, x2, x3p) = h[D3Q27System::INV_BW];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_BE, x1, x2, x3p) = h[D3Q27System::INV_BE];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_BS, x1, x2p, x3p) = h[D3Q27System::INV_BS];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_BN, x1, x2, x3p) = h[D3Q27System::INV_BN];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_BSW, x1p, x2p, x3p) = h[D3Q27System::INV_BSW];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_BSE, x1, x2p, x3p) = h[D3Q27System::INV_BSE];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_BNW, x1p, x2, x3p) = h[D3Q27System::INV_BNW];
+ // (*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p) = h[D3Q27System::INV_BNE];
+
+ // (*this->zeroDistributionsH)(x1, x2, x3) = h[D3Q27System::REST];
+
+ /////////////////// END OF OLD BGK SOLVER ///////////////////////////////
}
}
}
--
GitLab