diff --git a/src/gpu/VirtualFluids_GPU/GPU/Cumulant27chim.cu b/src/gpu/VirtualFluids_GPU/GPU/Cumulant27chim.cu
index a8c4b7ce64c98974cee06fef45f154f0eb179b73..49e2c10a2072762fc5568d6875ec4b9e8dcb7567 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/Cumulant27chim.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/Cumulant27chim.cu
@@ -107,27 +107,30 @@ __device__ Distributions27 getDistributions27(real* distributions, unsigned int
return dist;
}
-struct ABC
+struct DistributionWrapper
{
- __device__ ABC(uint k,
- uint kw,
- uint ks,
- uint kb,
- uint ksw,
- uint kbw,
- uint kbs,
- uint kbsw) :
+ __device__ DistributionWrapper(
+ real* distributions,
+ unsigned int size_Mat,
+ bool isEvenTimestep,
+ uint k,
+ uint* neighborX,
+ uint* neighborY,
+ uint* neighborZ) :
+ dist(getDistributions27(distributions, size_Mat, isEvenTimestep)),
k(k),
- kw (kw),
- ks (ks),
- kb (kb),
- ksw (ksw),
- kbw (kbw),
- kbs (kbs),
- kbsw (kbsw)
- {}
-
- __device__ void read(const Distributions27& dist)
+ kw (neighborX[k]),
+ ks (neighborY[k]),
+ kb (neighborZ[k]),
+ ksw (neighborY[kw]),
+ kbw (neighborZ[kw]),
+ kbs (neighborZ[ks]),
+ kbsw(neighborZ[ksw])
+ {
+ read();
+ }
+
+ __device__ void read()
{
distribution.f[VF::LBM::DIR::PZZ] = (dist.f[dirE ])[k];
distribution.f[VF::LBM::DIR::MZZ] = (dist.f[dirW ])[kw];
@@ -152,13 +155,13 @@ struct ABC
distribution.f[VF::LBM::DIR::PMP] = (dist.f[dirTSE ])[ks];
distribution.f[VF::LBM::DIR::MMP] = (dist.f[dirTSW ])[ksw];
distribution.f[VF::LBM::DIR::PPM] = (dist.f[dirBNE ])[kb];
- distribution.f[VF::LBM::DIR::MPM] = (dist.f[dirBNE ])[kb];
+ distribution.f[VF::LBM::DIR::MPM] = (dist.f[dirBNW ])[kbw];
distribution.f[VF::LBM::DIR::PMM] = (dist.f[dirBSE ])[kbs];
distribution.f[VF::LBM::DIR::MMM] = (dist.f[dirBSW ])[kbsw];
distribution.f[VF::LBM::DIR::ZZZ] = (dist.f[dirREST])[k];
}
- __device__ void write(Distributions27& dist)
+ __device__ void write()
{
(dist.f[dirE ])[k] = distribution.f[VF::LBM::DIR::PZZ];
(dist.f[dirW ])[kw] = distribution.f[VF::LBM::DIR::MZZ];
@@ -183,12 +186,14 @@ struct ABC
(dist.f[dirTSE ])[ks] = distribution.f[VF::LBM::DIR::PMP];
(dist.f[dirTSW ])[ksw] = distribution.f[VF::LBM::DIR::MMP];
(dist.f[dirBNE ])[kb] = distribution.f[VF::LBM::DIR::PPM];
- (dist.f[dirBNE ])[kb] = distribution.f[VF::LBM::DIR::MPM];
+ (dist.f[dirBNW ])[kbw] = distribution.f[VF::LBM::DIR::MPM];
(dist.f[dirBSE ])[kbs] = distribution.f[VF::LBM::DIR::PMM];
(dist.f[dirBSW ])[kbsw] = distribution.f[VF::LBM::DIR::MMM];
(dist.f[dirREST])[k] = distribution.f[VF::LBM::DIR::ZZZ];
}
+ Distributions27 dist;
+
VF::LBM::Distribution27 distribution;
const uint k;
@@ -201,7 +206,39 @@ struct ABC
const uint kbsw;
};
-////////////////////////////////////////////////////////////////////////////////
+__device__ unsigned int getNodeIndex()
+{
+ const unsigned x = threadIdx.x;
+ const unsigned y = blockIdx.x;
+ const unsigned z = blockIdx.y;
+
+ const unsigned nx = blockDim.x;
+ const unsigned ny = gridDim.x;
+
+ return nx*(ny*z + y) + x;
+}
+
+__device__ bool isValidFluidNode(uint k, int size_Mat, uint nodeType)
+{
+ return (k < size_Mat) && (nodeType == GEO_FLUID);
+}
+
+__device__ void getLevelForce(real fx, real fy, real fz, int level, real* forces)
+{
+ real fx_t {1.}, fy_t {1.}, fz_t {1.};
+ for (int i = 0; i < level; i++)
+ {
+ fx_t *= c2o1;
+ fy_t *= c2o1;
+ fz_t *= c2o1;
+ }
+
+ forces[0] = fx / fx_t;
+ forces[1] = fy / fy_t;
+ forces[2] = fz / fz_t;
+}
+
+
extern "C" __global__ void Cumulant_K17_LBM_Device_Kernel(
real omega,
uint* typeOfGridNode,
@@ -213,6 +250,39 @@ extern "C" __global__ void Cumulant_K17_LBM_Device_Kernel(
int level,
real* forces,
bool isEvenTimestep)
+{
+ const uint k = getNodeIndex();
+ const uint nodeType = typeOfGridNode[k];
+
+ if (isValidFluidNode(k, size_Mat, nodeType))
+ {
+ DistributionWrapper distributionWrapper {
+ distributions, size_Mat, isEvenTimestep, k, neighborX, neighborY, neighborZ
+ };
+
+ real level_forces[3];
+ getLevelForce(forces[0], forces[1], forces[2], level, level_forces);
+
+ VF::LBM::cumulantChimera(distributionWrapper.distribution, omega, level_forces);
+
+ distributionWrapper.write();
+ }
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+extern "C" __global__ void Cumulant_K17_LBM_Device_Kernel_old(
+ real omega,
+ uint* typeOfGridNode,
+ uint* neighborX,
+ uint* neighborY,
+ uint* neighborZ,
+ real* distributions,
+ int size_Mat,
+ int level,
+ real* forces,
+ bool isEvenTimestep)
{
//////////////////////////////////////////////////////////////////////////
//! Cumulant K17 Kernel is based on \ref
@@ -244,112 +314,555 @@ extern "C" __global__ void Cumulant_K17_LBM_Device_Kernel(
//!
Distributions27 dist = getDistributions27(distributions, size_Mat, isEvenTimestep);
- // ////////////////////////////////////////////////////////////////////////////////
- // //! - Set neighbor indices (necessary for indirect addressing)
- const uint kw = neighborX[k];
- const uint ks = neighborY[k];
- const uint kb = neighborZ[k];
- const uint ksw = neighborY[kw];
- const uint kbw = neighborZ[kw];
- const uint kbs = neighborZ[ks];
- const uint kbsw = neighborZ[ksw];
-
- ABC abc = {
- k,
- kw ,
- ks ,
- kb ,
- ksw ,
- kbw ,
- kbs ,
- kbsw
- };
+ ////////////////////////////////////////////////////////////////////////////////
+ //! - Set neighbor indices (necessary for indirect addressing)
+ uint kw = neighborX[k];
+ uint ks = neighborY[k];
+ uint kb = neighborZ[k];
+ uint ksw = neighborY[kw];
+ uint kbw = neighborZ[kw];
+ uint kbs = neighborZ[ks];
+ uint kbsw = neighborZ[ksw];
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Set local distributions
+ //!
+ // real mfcbb = (dist.f[dirE ])[k];
+ // real mfabb = (dist.f[dirW ])[kw];
+ // real mfbcb = (dist.f[dirN ])[k];
+ // real mfbab = (dist.f[dirS ])[ks];
+ // real mfbbc = (dist.f[dirT ])[k];
+ // real mfbba = (dist.f[dirB ])[kb];
+ // real mfccb = (dist.f[dirNE ])[k];
+ // real mfaab = (dist.f[dirSW ])[ksw];
+ // real mfcab = (dist.f[dirSE ])[ks];
+ // real mfacb = (dist.f[dirNW ])[kw];
+ // real mfcbc = (dist.f[dirTE ])[k];
+ // real mfaba = (dist.f[dirBW ])[kbw];
+ // real mfcba = (dist.f[dirBE ])[kb];
+ // real mfabc = (dist.f[dirTW ])[kw];
+ // real mfbcc = (dist.f[dirTN ])[k];
+ // real mfbaa = (dist.f[dirBS ])[kbs];
+ // real mfbca = (dist.f[dirBN ])[kb];
+ // real mfbac = (dist.f[dirTS ])[ks];
+ // real mfbbb = (dist.f[dirREST])[k];
+ // real mfccc = (dist.f[dirTNE ])[k];
+ // real mfaac = (dist.f[dirTSW ])[ksw];
+ // real mfcac = (dist.f[dirTSE ])[ks];
+ // real mfacc = (dist.f[dirTNW ])[kw];
+ // real mfcca = (dist.f[dirBNE ])[kb];
+ // real mfaaa = (dist.f[dirBSW ])[kbsw];
+ // real mfcaa = (dist.f[dirBSE ])[kbs];
+ // real mfaca = (dist.f[dirBNW ])[kbw];
+
+ VF::LBM::Distribution27 distribution_in;
+
+ distribution_in.f[VF::LBM::DIR::PZZ] = (dist.f[dirE ])[k];
+ distribution_in.f[VF::LBM::DIR::MZZ] = (dist.f[dirW ])[kw];
+ distribution_in.f[VF::LBM::DIR::ZPZ] = (dist.f[dirN ])[k];
+ distribution_in.f[VF::LBM::DIR::ZMZ] = (dist.f[dirS ])[ks];
+ distribution_in.f[VF::LBM::DIR::ZZP] = (dist.f[dirT ])[k];
+ distribution_in.f[VF::LBM::DIR::ZZM] = (dist.f[dirB ])[kb];
+ distribution_in.f[VF::LBM::DIR::PPZ] = (dist.f[dirNE ])[k];
+ distribution_in.f[VF::LBM::DIR::MMZ] = (dist.f[dirSW ])[ksw];
+ distribution_in.f[VF::LBM::DIR::PMZ] = (dist.f[dirSE ])[ks];
+ distribution_in.f[VF::LBM::DIR::MPZ] = (dist.f[dirNW ])[kw];
+ distribution_in.f[VF::LBM::DIR::PZP] = (dist.f[dirTE ])[k];
+ distribution_in.f[VF::LBM::DIR::MZM] = (dist.f[dirBW ])[kbw];
+ distribution_in.f[VF::LBM::DIR::PZM] = (dist.f[dirBE ])[kb];
+ distribution_in.f[VF::LBM::DIR::MZP] = (dist.f[dirTW ])[kw];
+ distribution_in.f[VF::LBM::DIR::ZPP] = (dist.f[dirTN ])[k];
+ distribution_in.f[VF::LBM::DIR::ZMM] = (dist.f[dirBS ])[kbs];
+ distribution_in.f[VF::LBM::DIR::ZPM] = (dist.f[dirBN ])[kb];
+ distribution_in.f[VF::LBM::DIR::ZMP] = (dist.f[dirTS ])[ks];
+ distribution_in.f[VF::LBM::DIR::PPP] = (dist.f[dirTNE ])[k];
+ distribution_in.f[VF::LBM::DIR::MPP] = (dist.f[dirTNW ])[kw];
+ distribution_in.f[VF::LBM::DIR::PMP] = (dist.f[dirTSE ])[ks];
+ distribution_in.f[VF::LBM::DIR::MMP] = (dist.f[dirTSW ])[ksw];
+ distribution_in.f[VF::LBM::DIR::PPM] = (dist.f[dirBNE ])[kb];
+ distribution_in.f[VF::LBM::DIR::MPM] = (dist.f[dirBNW ])[kbw];
+ distribution_in.f[VF::LBM::DIR::PMM] = (dist.f[dirBSE ])[kbs];
+ distribution_in.f[VF::LBM::DIR::MMM] = (dist.f[dirBSW ])[kbsw];
+ distribution_in.f[VF::LBM::DIR::ZZZ] = (dist.f[dirREST])[k];
+
+ real mfcbb = distribution_in.f[VF::LBM::DIR::PZZ];
+ real mfabb = distribution_in.f[VF::LBM::DIR::MZZ];
+ real mfbcb = distribution_in.f[VF::LBM::DIR::ZPZ];
+ real mfbab = distribution_in.f[VF::LBM::DIR::ZMZ];
+ real mfbbc = distribution_in.f[VF::LBM::DIR::ZZP];
+ real mfbba = distribution_in.f[VF::LBM::DIR::ZZM];
+ real mfccb = distribution_in.f[VF::LBM::DIR::PPZ];
+ real mfaab = distribution_in.f[VF::LBM::DIR::MMZ];
+ real mfcab = distribution_in.f[VF::LBM::DIR::PMZ];
+ real mfacb = distribution_in.f[VF::LBM::DIR::MPZ];
+ real mfcbc = distribution_in.f[VF::LBM::DIR::PZP];
+ real mfaba = distribution_in.f[VF::LBM::DIR::MZM];
+ real mfcba = distribution_in.f[VF::LBM::DIR::PZM];
+ real mfabc = distribution_in.f[VF::LBM::DIR::MZP];
+ real mfbcc = distribution_in.f[VF::LBM::DIR::ZPP];
+ real mfbaa = distribution_in.f[VF::LBM::DIR::ZMM];
+ real mfbca = distribution_in.f[VF::LBM::DIR::ZPM];
+ real mfbac = distribution_in.f[VF::LBM::DIR::ZMP];
+ real mfccc = distribution_in.f[VF::LBM::DIR::PPP];
+ real mfacc = distribution_in.f[VF::LBM::DIR::MPP];
+ real mfcac = distribution_in.f[VF::LBM::DIR::PMP];
+ real mfaac = distribution_in.f[VF::LBM::DIR::MMP];
+ real mfcca = distribution_in.f[VF::LBM::DIR::PPM];
+ real mfaca = distribution_in.f[VF::LBM::DIR::MPM];
+ real mfcaa = distribution_in.f[VF::LBM::DIR::PMM];
+ real mfaaa = distribution_in.f[VF::LBM::DIR::MMM];
+ real mfbbb = distribution_in.f[VF::LBM::DIR::ZZZ];
+
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Calculate density and velocity using pyramid summation for low round-off errors as in Eq. (J1)-(J3) \ref
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
+ //!
+ real drho =
+ ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) +
+ (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
+ ((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb;
+ real rho = c1o1 + drho;
+ real OOrho = c1o1 / rho;
+ real vvx =
+ ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+ (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
+ (mfcbb - mfabb)) * OOrho;
+ real vvy =
+ ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+ (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
+ (mfbcb - mfbab)) * OOrho;
+ real vvz =
+ ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+ (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
+ (mfbbc - mfbba)) * OOrho;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Add half of the acceleration (body force) to the velocity as in Eq. (42) \ref
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
+ //!
+ real fx = forces[0];
+ real fy = forces[1];
+ real fz = forces[2];
+
+ real fx_t {1.}, fy_t {1.}, fz_t {1.};
+ for (int i = 0; i < level; i++)
+ {
+ fx_t *= c2o1;
+ fy_t *= c2o1;
+ fz_t *= c2o1;
+ }
+
+ fx /= fx_t;
+ fy /= fy_t;
+ fz /= fz_t;
+ //real forces[3] {fx, fy, fz};
- // ////////////////////////////////////////////////////////////////////////////////////
- // //! - Set local distributions
- // //!
- // VF::LBM::Distribution27 distribution_in;
-
- // distribution_in.f[VF::LBM::DIR::PZZ] = (dist.f[dirE ])[k];
- // distribution_in.f[VF::LBM::DIR::MZZ] = (dist.f[dirW ])[kw];
- // distribution_in.f[VF::LBM::DIR::ZPZ] = (dist.f[dirN ])[k];
- // distribution_in.f[VF::LBM::DIR::ZMZ] = (dist.f[dirS ])[ks];
- // distribution_in.f[VF::LBM::DIR::ZZP] = (dist.f[dirT ])[k];
- // distribution_in.f[VF::LBM::DIR::ZZM] = (dist.f[dirB ])[kb];
- // distribution_in.f[VF::LBM::DIR::PPZ] = (dist.f[dirNE ])[k];
- // distribution_in.f[VF::LBM::DIR::MMZ] = (dist.f[dirSW ])[ksw];
- // distribution_in.f[VF::LBM::DIR::PMZ] = (dist.f[dirSE ])[ks];
- // distribution_in.f[VF::LBM::DIR::MPZ] = (dist.f[dirNW ])[kw];
- // distribution_in.f[VF::LBM::DIR::PZP] = (dist.f[dirTE ])[k];
- // distribution_in.f[VF::LBM::DIR::MZM] = (dist.f[dirBW ])[kbw];
- // distribution_in.f[VF::LBM::DIR::PZM] = (dist.f[dirBE ])[kb];
- // distribution_in.f[VF::LBM::DIR::MZP] = (dist.f[dirTW ])[kw];
- // distribution_in.f[VF::LBM::DIR::ZPP] = (dist.f[dirTN ])[k];
- // distribution_in.f[VF::LBM::DIR::ZMM] = (dist.f[dirBS ])[kbs];
- // distribution_in.f[VF::LBM::DIR::ZPM] = (dist.f[dirBN ])[kb];
- // distribution_in.f[VF::LBM::DIR::ZMP] = (dist.f[dirTS ])[ks];
- // distribution_in.f[VF::LBM::DIR::PPP] = (dist.f[dirTNE ])[k];
- // distribution_in.f[VF::LBM::DIR::MPP] = (dist.f[dirTNW ])[kw];
- // distribution_in.f[VF::LBM::DIR::PMP] = (dist.f[dirTSE ])[ks];
- // distribution_in.f[VF::LBM::DIR::MMP] = (dist.f[dirTSW ])[ksw];
- // distribution_in.f[VF::LBM::DIR::PPM] = (dist.f[dirBNE ])[kb];
- // distribution_in.f[VF::LBM::DIR::MPM] = (dist.f[dirBNE ])[kb];
- // distribution_in.f[VF::LBM::DIR::PMM] = (dist.f[dirBSE ])[kbs];
- // distribution_in.f[VF::LBM::DIR::MMM] = (dist.f[dirBSW ])[kbsw];
- // distribution_in.f[VF::LBM::DIR::ZZZ] = (dist.f[dirREST])[k];
-
- abc.read(dist);
-
- real vx = forces[0] / pow((double)c2o1, (double)level);
- real vy = forces[1] / pow((double)c2o1, (double)level);
- real vz = forces[2] / pow((double)c2o1, (double)level);
- real forces[3] {vx, vy, vz};
-
- // VF::LBM::Distribution27 distribution_out = VF::LBM::cumulantChimera(distribution_in, omega, forces);
-
- abc.distribution = VF::LBM::cumulantChimera(abc.distribution, omega, forces);
-
- abc.write(dist);
+ vvx += fx * c1o2;
+ vvy += fy * c1o2;
+ vvz += fz * c1o2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ // calculate the square of velocities for this lattice node
+ real vx2 = vvx*vvx;
+ real vy2 = vvy*vvy;
+ real vz2 = vvz*vvz;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Set relaxation limiters for third order cumulants to default value \f$ \lambda=0.001 \f$ according to section 6 in \ref
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //!
+ real wadjust;
+ real qudricLimitP = c1o100;
+ real qudricLimitM = c1o100;
+ real qudricLimitD = c1o100;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Chimera transform from well conditioned distributions to central moments as defined in Appendix J in \ref
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
+ //! see also Eq. (6)-(14) in \ref
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ VF::LBM::forwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36o1, c1o36);
+ VF::LBM::forwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9o1, c1o9);
+ VF::LBM::forwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36o1, c1o36);
+ VF::LBM::forwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9o1, c1o9);
+ VF::LBM::forwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
+ VF::LBM::forwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9o1, c1o9);
+ VF::LBM::forwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36o1, c1o36);
+ VF::LBM::forwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9o1, c1o9);
+ VF::LBM::forwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36o1, c1o36);
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ VF::LBM::forwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6o1, c1o6);
+ VF::LBM::forwardChimera( mfaab, mfabb, mfacb, vvy, vy2);
+ VF::LBM::forwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18o1, c1o18);
+ VF::LBM::forwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
+ VF::LBM::forwardChimera( mfbab, mfbbb, mfbcb, vvy, vy2);
+ VF::LBM::forwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
+ VF::LBM::forwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6o1, c1o6);
+ VF::LBM::forwardChimera( mfcab, mfcbb, mfccb, vvy, vy2);
+ VF::LBM::forwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18o1, c1o18);
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ VF::LBM::forwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1o1, c1o1);
+ VF::LBM::forwardChimera( mfaba, mfbba, mfcba, vvx, vx2);
+ VF::LBM::forwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3o1, c1o3);
+ VF::LBM::forwardChimera( mfaab, mfbab, mfcab, vvx, vx2);
+ VF::LBM::forwardChimera( mfabb, mfbbb, mfcbb, vvx, vx2);
+ VF::LBM::forwardChimera( mfacb, mfbcb, mfccb, vvx, vx2);
+ VF::LBM::forwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3o1, c1o3);
+ VF::LBM::forwardChimera( mfabc, mfbbc, mfcbc, vvx, vx2);
+ VF::LBM::forwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c3o1, c1o9);
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Setting relaxation rates for non-hydrodynamic cumulants (default values). Variable names and equations according to
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //! => [NAME IN PAPER]=[NAME IN CODE]=[DEFAULT VALUE].
+ //! - Trace of second order cumulants \f$ C_{200}+C_{020}+C_{002} \f$ used to adjust bulk viscosity:\f$\omega_2=OxxPyyPzz=1.0 \f$.
+ //! - Third order cumulants \f$ C_{120}+C_{102}, C_{210}+C_{012}, C_{201}+C_{021} \f$: \f$ \omega_3=OxyyPxzz \f$ set according to Eq. (111) with simplifications assuming \f$ \omega_2=1.0\f$.
+ //! - Third order cumulants \f$ C_{120}-C_{102}, C_{210}-C_{012}, C_{201}-C_{021} \f$: \f$ \omega_4 = OxyyMxzz \f$ set according to Eq. (112) with simplifications assuming \f$ \omega_2 = 1.0\f$.
+ //! - Third order cumulants \f$ C_{111} \f$: \f$ \omega_5 = Oxyz \f$ set according to Eq. (113) with simplifications assuming \f$ \omega_2 = 1.0\f$ (modify for different bulk viscosity).
+ //! - Fourth order cumulants \f$ C_{220}, C_{202}, C_{022}, C_{211}, C_{121}, C_{112} \f$: for simplification all set to the same default value \f$ \omega_6=\omega_7=\omega_8=O4=1.0 \f$.
+ //! - Fifth order cumulants \f$ C_{221}, C_{212}, C_{122}\f$: \f$\omega_9=O5=1.0\f$.
+ //! - Sixth order cumulant \f$ C_{222}\f$: \f$\omega_{10}=O6=1.0\f$.
+ //!
+ ////////////////////////////////////////////////////////////
+ //2.
+ real OxxPyyPzz = c1o1;
+ ////////////////////////////////////////////////////////////
+ //3.
+ real OxyyPxzz = c8o1 * (-c2o1 + omega) * ( c1o1 + c2o1*omega) / (-c8o1 - c14o1*omega + c7o1*omega*omega);
+ real OxyyMxzz = c8o1 * (-c2o1 + omega) * (-c7o1 + c4o1*omega) / (c56o1 - c50o1*omega + c9o1*omega*omega);
+ real Oxyz = c24o1 * (-c2o1 + omega) * (-c2o1 - c7o1*omega + c3o1*omega*omega) / (c48o1 + c152o1*omega - c130o1*omega*omega + c29o1*omega*omega*omega);
+ ////////////////////////////////////////////////////////////
+ //4.
+ real O4 = c1o1;
+ ////////////////////////////////////////////////////////////
+ //5.
+ real O5 = c1o1;
+ ////////////////////////////////////////////////////////////
+ //6.
+ real O6 = c1o1;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - A and B: parameters for fourth order convergence of the diffusion term according to Eq. (114) and (115)
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //! with simplifications assuming \f$ \omega_2 = 1.0 \f$ (modify for different bulk viscosity).
+ //!
+ real A = (c4o1 + c2o1*omega - c3o1*omega*omega) / (c2o1 - c7o1*omega + c5o1*omega*omega);
+ real B = (c4o1 + c28o1*omega - c14o1*omega*omega) / (c6o1 - c21o1*omega + c15o1*omega*omega);
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Compute cumulants from central moments according to Eq. (20)-(23) in
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////
+ //4.
+ real CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + c2o1 * mfbba * mfbab) * OOrho;
+ real CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + c2o1 * mfbba * mfabb) * OOrho;
+ real CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + c2o1 * mfbab * mfabb) * OOrho;
+ real CUMcca = mfcca - (((mfcaa * mfaca + c2o1 * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) * OOrho - c1o9*(drho * OOrho));
+ real CUMcac = mfcac - (((mfcaa * mfaac + c2o1 * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) * OOrho - c1o9*(drho * OOrho));
+ real CUMacc = mfacc - (((mfaac * mfaca + c2o1 * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) * OOrho - c1o9*(drho * OOrho));
+ ////////////////////////////////////////////////////////////
+ //5.
+ real CUMbcc = mfbcc - ((mfaac * mfbca + mfaca * mfbac + c4o1 * mfabb * mfbbb + c2o1 * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac)) * OOrho;
+ real CUMcbc = mfcbc - ((mfaac * mfcba + mfcaa * mfabc + c4o1 * mfbab * mfbbb + c2o1 * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc)) * OOrho;
+ real CUMccb = mfccb - ((mfcaa * mfacb + mfaca * mfcab + c4o1 * mfbba * mfbbb + c2o1 * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab)) * OOrho;
+ ////////////////////////////////////////////////////////////
+ //6.
+ real CUMccc = mfccc + ((-c4o1 * mfbbb * mfbbb
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - c4o1 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+ - c2o1 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
+ + (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + c2o1 * (mfcaa * mfaca * mfaac)
+ + c16o1 * mfbba * mfbab * mfabb) * OOrho * OOrho
+ - c1o3 * (mfacc + mfcac + mfcca) * OOrho
+ - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
+ + (c2o1 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
+ + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) * OOrho * OOrho * c2o3
+ + c1o27*((drho * drho - drho) * OOrho * OOrho));
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Compute linear combinations of second and third order cumulants
+ //!
+ ////////////////////////////////////////////////////////////
+ //2.
+ real mxxPyyPzz = mfcaa + mfaca + mfaac;
+ real mxxMyy = mfcaa - mfaca;
+ real mxxMzz = mfcaa - mfaac;
+ ////////////////////////////////////////////////////////////
+ //3.
+ real mxxyPyzz = mfcba + mfabc;
+ real mxxyMyzz = mfcba - mfabc;
+ real mxxzPyyz = mfcab + mfacb;
+ real mxxzMyyz = mfcab - mfacb;
+ real mxyyPxzz = mfbca + mfbac;
+ real mxyyMxzz = mfbca - mfbac;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //incl. correction
+ ////////////////////////////////////////////////////////////
+ //! - Compute velocity gradients from second order cumulants according to Eq. (27)-(32)
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //! Further explanations of the correction in viscosity in Appendix H of
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
+ //! Note that the division by rho is omitted here as we need rho times the gradients later.
+ //!
+ real Dxy = -c3o1*omega*mfbba;
+ real Dxz = -c3o1*omega*mfbab;
+ real Dyz = -c3o1*omega*mfabb;
+ real dxux = c1o2 * (-omega) *(mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
+ real dyuy = dxux + omega * c3o2 * mxxMyy;
+ real dzuz = dxux + omega * c3o2 * mxxMzz;
+ ////////////////////////////////////////////////////////////
+ //! - Relaxation of second order cumulants with correction terms according to Eq. (33)-(35) in
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //!
+ mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz) - c3o1 * (c1o1 - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
+ mxxMyy += omega * (-mxxMyy) - c3o1 * (c1o1 + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
+ mxxMzz += omega * (-mxxMzz) - c3o1 * (c1o1 + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////no correction
+ //mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz);
+ //mxxMyy += -(-omega) * (-mxxMyy);
+ //mxxMzz += -(-omega) * (-mxxMzz);
+ //////////////////////////////////////////////////////////////////////////
+ mfabb += omega * (-mfabb);
+ mfbab += omega * (-mfbab);
+ mfbba += omega * (-mfbba);
+ ////////////////////////////////////////////////////////////////////////////////////
+ //relax
+ //////////////////////////////////////////////////////////////////////////
+ // incl. limiter
+ //! - Relaxation of third order cumulants including limiter according to Eq. (116)-(123)
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //!
+ wadjust = Oxyz + (c1o1 - Oxyz)*abs(mfbbb) / (abs(mfbbb) + qudricLimitD);
+ mfbbb += wadjust * (-mfbbb);
+ wadjust = OxyyPxzz + (c1o1 - OxyyPxzz)*abs(mxxyPyzz) / (abs(mxxyPyzz) + qudricLimitP);
+ mxxyPyzz += wadjust * (-mxxyPyzz);
+ wadjust = OxyyMxzz + (c1o1 - OxyyMxzz)*abs(mxxyMyzz) / (abs(mxxyMyzz) + qudricLimitM);
+ mxxyMyzz += wadjust * (-mxxyMyzz);
+ wadjust = OxyyPxzz + (c1o1 - OxyyPxzz)*abs(mxxzPyyz) / (abs(mxxzPyyz) + qudricLimitP);
+ mxxzPyyz += wadjust * (-mxxzPyyz);
+ wadjust = OxyyMxzz + (c1o1 - OxyyMxzz)*abs(mxxzMyyz) / (abs(mxxzMyyz) + qudricLimitM);
+ mxxzMyyz += wadjust * (-mxxzMyyz);
+ wadjust = OxyyPxzz + (c1o1 - OxyyPxzz)*abs(mxyyPxzz) / (abs(mxyyPxzz) + qudricLimitP);
+ mxyyPxzz += wadjust * (-mxyyPxzz);
+ wadjust = OxyyMxzz + (c1o1 - OxyyMxzz)*abs(mxyyMxzz) / (abs(mxyyMxzz) + qudricLimitM);
+ mxyyMxzz += wadjust * (-mxyyMxzz);
+ //////////////////////////////////////////////////////////////////////////
+ // no limiter
+ //mfbbb += OxyyMxzz * (-mfbbb);
+ //mxxyPyzz += OxyyPxzz * (-mxxyPyzz);
+ //mxxyMyzz += OxyyMxzz * (-mxxyMyzz);
+ //mxxzPyyz += OxyyPxzz * (-mxxzPyyz);
+ //mxxzMyyz += OxyyMxzz * (-mxxzMyyz);
+ //mxyyPxzz += OxyyPxzz * (-mxyyPxzz);
+ //mxyyMxzz += OxyyMxzz * (-mxyyMxzz);
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Compute inverse linear combinations of second and third order cumulants
+ //!
+ mfcaa = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz);
+ mfaca = c1o3 * (-c2o1* mxxMyy + mxxMzz + mxxPyyPzz);
+ mfaac = c1o3 * (mxxMyy - c2o1* mxxMzz + mxxPyyPzz);
+ mfcba = ( mxxyMyzz + mxxyPyzz) * c1o2;
+ mfabc = (-mxxyMyzz + mxxyPyzz) * c1o2;
+ mfcab = ( mxxzMyyz + mxxzPyyz) * c1o2;
+ mfacb = (-mxxzMyyz + mxxzPyyz) * c1o2;
+ mfbca = ( mxyyMxzz + mxyyPxzz) * c1o2;
+ mfbac = (-mxyyMxzz + mxyyPxzz) * c1o2;
+ //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ //4.
+ // no limiter
+ //! - Relax fourth order cumulants to modified equilibrium for fourth order convergence of diffusion according to Eq. (43)-(48)
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //!
+ CUMacc = -O4*(c1o1 / omega - c1o2) * (dyuy + dzuz) * c2o3 * A + (c1o1 - O4) * (CUMacc);
+ CUMcac = -O4*(c1o1 / omega - c1o2) * (dxux + dzuz) * c2o3 * A + (c1o1 - O4) * (CUMcac);
+ CUMcca = -O4*(c1o1 / omega - c1o2) * (dyuy + dxux) * c2o3 * A + (c1o1 - O4) * (CUMcca);
+ CUMbbc = -O4*(c1o1 / omega - c1o2) * Dxy * c1o3 * B + (c1o1 - O4) * (CUMbbc);
+ CUMbcb = -O4*(c1o1 / omega - c1o2) * Dxz * c1o3 * B + (c1o1 - O4) * (CUMbcb);
+ CUMcbb = -O4*(c1o1 / omega - c1o2) * Dyz * c1o3 * B + (c1o1 - O4) * (CUMcbb);
+ //////////////////////////////////////////////////////////////////////////
+ //5.
+ CUMbcc += O5 * (-CUMbcc);
+ CUMcbc += O5 * (-CUMcbc);
+ CUMccb += O5 * (-CUMccb);
+ //////////////////////////////////////////////////////////////////////////
+ //6.
+ CUMccc += O6 * (-CUMccc);
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Compute central moments from post collision cumulants according to Eq. (53)-(56) in
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //!
+ //////////////////////////////////////////////////////////////////////////
+ //4.
+ mfcbb = CUMcbb + c1o3*((c3o1*mfcaa + c1o1) * mfabb + c6o1 * mfbba * mfbab) * OOrho;
+ mfbcb = CUMbcb + c1o3*((c3o1*mfaca + c1o1) * mfbab + c6o1 * mfbba * mfabb) * OOrho;
+ mfbbc = CUMbbc + c1o3*((c3o1*mfaac + c1o1) * mfbba + c6o1 * mfbab * mfabb) * OOrho;
+ mfcca = CUMcca + (((mfcaa * mfaca + c2o1 * mfbba * mfbba)*c9o1 + c3o1 * (mfcaa + mfaca)) * OOrho - (drho * OOrho))*c1o9;
+ mfcac = CUMcac + (((mfcaa * mfaac + c2o1 * mfbab * mfbab)*c9o1 + c3o1 * (mfcaa + mfaac)) * OOrho - (drho * OOrho))*c1o9;
+ mfacc = CUMacc + (((mfaac * mfaca + c2o1 * mfabb * mfabb)*c9o1 + c3o1 * (mfaac + mfaca)) * OOrho - (drho * OOrho))*c1o9;
+ //////////////////////////////////////////////////////////////////////////
+ //5.
+ mfbcc = CUMbcc + c1o3 *(c3o1*(mfaac * mfbca + mfaca * mfbac + c4o1 * mfabb * mfbbb + c2o1 * (mfbab * mfacb + mfbba * mfabc)) + (mfbca + mfbac)) * OOrho;
+ mfcbc = CUMcbc + c1o3 *(c3o1*(mfaac * mfcba + mfcaa * mfabc + c4o1 * mfbab * mfbbb + c2o1 * (mfabb * mfcab + mfbba * mfbac)) + (mfcba + mfabc)) * OOrho;
+ mfccb = CUMccb + c1o3 *(c3o1*(mfcaa * mfacb + mfaca * mfcab + c4o1 * mfbba * mfbbb + c2o1 * (mfbab * mfbca + mfabb * mfcba)) + (mfacb + mfcab)) * OOrho;
+ //////////////////////////////////////////////////////////////////////////
+ //6.
+ mfccc = CUMccc - ((-c4o1 * mfbbb * mfbbb
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - c4o1 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+ - c2o1 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
+ + (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + c2o1 * (mfcaa * mfaca * mfaac)
+ + c16o1 * mfbba * mfbab * mfabb) * OOrho * OOrho
+ - c1o3 * (mfacc + mfcac + mfcca) * OOrho
+ - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
+ + (c2o1 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
+ + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) * OOrho * OOrho * c2o3
+ + c1o27*((drho * drho - drho) * OOrho * OOrho));
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Add acceleration (body force) to first order cumulants according to Eq. (85)-(87) in
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
+ //!
+ mfbaa = -mfbaa;
+ mfaba = -mfaba;
+ mfaab = -mfaab;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Chimera transform from central moments to well conditioned distributions as defined in Appendix J in
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
+ //! see also Eq. (88)-(96) in
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ VF::LBM::backwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1o1, c1o1);
+ VF::LBM::backwardChimera( mfaba, mfbba, mfcba, vvx, vx2);
+ VF::LBM::backwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3o1, c1o3);
+ VF::LBM::backwardChimera( mfaab, mfbab, mfcab, vvx, vx2);
+ VF::LBM::backwardChimera( mfabb, mfbbb, mfcbb, vvx, vx2);
+ VF::LBM::backwardChimera( mfacb, mfbcb, mfccb, vvx, vx2);
+ VF::LBM::backwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3o1, c1o3);
+ VF::LBM::backwardChimera( mfabc, mfbbc, mfcbc, vvx, vx2);
+ VF::LBM::backwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c9o1, c1o9);
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ VF::LBM::backwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6o1, c1o6);
+ VF::LBM::backwardChimera( mfaab, mfabb, mfacb, vvy, vy2);
+ VF::LBM::backwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18o1, c1o18);
+ VF::LBM::backwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
+ VF::LBM::backwardChimera( mfbab, mfbbb, mfbcb, vvy, vy2);
+ VF::LBM::backwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
+ VF::LBM::backwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6o1, c1o6);
+ VF::LBM::backwardChimera( mfcab, mfcbb, mfccb, vvy, vy2);
+ VF::LBM::backwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18o1, c1o18);
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ VF::LBM::backwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36o1, c1o36);
+ VF::LBM::backwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9o1, c1o9);
+ VF::LBM::backwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36o1, c1o36);
+ VF::LBM::backwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9o1, c1o9);
+ VF::LBM::backwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
+ VF::LBM::backwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9o1, c1o9);
+ VF::LBM::backwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36o1, c1o36);
+ VF::LBM::backwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9o1, c1o9);
+ VF::LBM::backwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36o1, c1o36);
////////////////////////////////////////////////////////////////////////////////////
//! - Write distributions: style of reading and writing the distributions from/to
//! stored arrays dependent on timestep is based on the esoteric twist algorithm
//! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
//!
- // (dist.f[dirE ])[k] = distribution_out.f[VF::LBM::DIR::PZZ];
- // (dist.f[dirW ])[kw] = distribution_out.f[VF::LBM::DIR::MZZ];
- // (dist.f[dirN ])[k] = distribution_out.f[VF::LBM::DIR::ZPZ];
- // (dist.f[dirS ])[ks] = distribution_out.f[VF::LBM::DIR::ZMZ];
- // (dist.f[dirT ])[k] = distribution_out.f[VF::LBM::DIR::ZZP];
- // (dist.f[dirB ])[kb] = distribution_out.f[VF::LBM::DIR::ZZM];
- // (dist.f[dirNE ])[k] = distribution_out.f[VF::LBM::DIR::PPZ];
- // (dist.f[dirSW ])[ksw] = distribution_out.f[VF::LBM::DIR::MMZ];
- // (dist.f[dirSE ])[ks] = distribution_out.f[VF::LBM::DIR::PMZ];
- // (dist.f[dirNW ])[kw] = distribution_out.f[VF::LBM::DIR::MPZ];
- // (dist.f[dirTE ])[k] = distribution_out.f[VF::LBM::DIR::PZP];
- // (dist.f[dirBW ])[kbw] = distribution_out.f[VF::LBM::DIR::MZM];
- // (dist.f[dirBE ])[kb] = distribution_out.f[VF::LBM::DIR::PZM];
- // (dist.f[dirTW ])[kw] = distribution_out.f[VF::LBM::DIR::MZP];
- // (dist.f[dirTN ])[k] = distribution_out.f[VF::LBM::DIR::ZPP];
- // (dist.f[dirBS ])[kbs] = distribution_out.f[VF::LBM::DIR::ZMM];
- // (dist.f[dirBN ])[kb] = distribution_out.f[VF::LBM::DIR::ZPM];
- // (dist.f[dirTS ])[ks] = distribution_out.f[VF::LBM::DIR::ZMP];
- // (dist.f[dirTNE ])[k] = distribution_out.f[VF::LBM::DIR::PPP];
- // (dist.f[dirTNW ])[kw] = distribution_out.f[VF::LBM::DIR::MPP];
- // (dist.f[dirTSE ])[ks] = distribution_out.f[VF::LBM::DIR::PMP];
- // (dist.f[dirTSW ])[ksw] = distribution_out.f[VF::LBM::DIR::MMP];
- // (dist.f[dirBNE ])[kb] = distribution_out.f[VF::LBM::DIR::PPM];
- // (dist.f[dirBNE ])[kb] = distribution_out.f[VF::LBM::DIR::MPM];
- // (dist.f[dirBSE ])[kbs] = distribution_out.f[VF::LBM::DIR::PMM];
- // (dist.f[dirBSW ])[kbsw] = distribution_out.f[VF::LBM::DIR::MMM];
- // (dist.f[dirREST])[k] = distribution_out.f[VF::LBM::DIR::ZZZ];
+ // (dist.f[dirE ])[k ] = mfabb;
+ // (dist.f[dirW ])[kw ] = mfcbb;
+ // (dist.f[dirN ])[k ] = mfbab;
+ // (dist.f[dirS ])[ks ] = mfbcb;
+ // (dist.f[dirT ])[k ] = mfbba;
+ // (dist.f[dirB ])[kb ] = mfbbc;
+ // (dist.f[dirNE ])[k ] = mfaab;
+ // (dist.f[dirSW ])[ksw ] = mfccb;
+ // (dist.f[dirSE ])[ks ] = mfacb;
+ // (dist.f[dirNW ])[kw ] = mfcab;
+ // (dist.f[dirTE ])[k ] = mfaba;
+ // (dist.f[dirBW ])[kbw ] = mfcbc;
+ // (dist.f[dirBE ])[kb ] = mfabc;
+ // (dist.f[dirTW ])[kw ] = mfcba;
+ // (dist.f[dirTN ])[k ] = mfbaa;
+ // (dist.f[dirBS ])[kbs ] = mfbcc;
+ // (dist.f[dirBN ])[kb ] = mfbac;
+ // (dist.f[dirTS ])[ks ] = mfbca;
+ // (dist.f[dirREST])[k ] = mfbbb;
+ // (dist.f[dirTNE ])[k ] = mfaaa;
+ // (dist.f[dirTSE ])[ks ] = mfaca;
+ // (dist.f[dirBNE ])[kb ] = mfaac;
+ // (dist.f[dirBSE ])[kbs ] = mfacc;
+ // (dist.f[dirTNW ])[kw ] = mfcaa;
+ // (dist.f[dirTSW ])[ksw ] = mfcca;
+ // (dist.f[dirBNW ])[kbw ] = mfcac;
+ // (dist.f[dirBSW ])[kbsw] = mfccc;
+
+ VF::LBM::Distribution27 distribution_out;
+
+ distribution_out.f[VF::LBM::DIR::MZZ] = mfcbb;
+ distribution_out.f[VF::LBM::DIR::PZZ] = mfabb;
+ distribution_out.f[VF::LBM::DIR::ZMZ] = mfbcb;
+ distribution_out.f[VF::LBM::DIR::ZPZ] = mfbab;
+ distribution_out.f[VF::LBM::DIR::ZZM] = mfbbc;
+ distribution_out.f[VF::LBM::DIR::ZZP] = mfbba;
+ distribution_out.f[VF::LBM::DIR::MMZ] = mfccb;
+ distribution_out.f[VF::LBM::DIR::PPZ] = mfaab;
+ distribution_out.f[VF::LBM::DIR::MPZ] = mfcab;
+ distribution_out.f[VF::LBM::DIR::PMZ] = mfacb;
+ distribution_out.f[VF::LBM::DIR::MZM] = mfcbc;
+ distribution_out.f[VF::LBM::DIR::PZP] = mfaba;
+ distribution_out.f[VF::LBM::DIR::MZP] = mfcba;
+ distribution_out.f[VF::LBM::DIR::PZM] = mfabc;
+ distribution_out.f[VF::LBM::DIR::ZMM] = mfbcc;
+ distribution_out.f[VF::LBM::DIR::ZPP] = mfbaa;
+ distribution_out.f[VF::LBM::DIR::ZMP] = mfbca;
+ distribution_out.f[VF::LBM::DIR::ZPM] = mfbac;
+ distribution_out.f[VF::LBM::DIR::MMM] = mfccc;
+ distribution_out.f[VF::LBM::DIR::PMM] = mfacc;
+ distribution_out.f[VF::LBM::DIR::MPM] = mfcac;
+ distribution_out.f[VF::LBM::DIR::PPM] = mfaac;
+ distribution_out.f[VF::LBM::DIR::MMP] = mfcca;
+ distribution_out.f[VF::LBM::DIR::PMP] = mfaca;
+ distribution_out.f[VF::LBM::DIR::MPP] = mfcaa;
+ distribution_out.f[VF::LBM::DIR::PPP] = mfaaa;
+ distribution_out.f[VF::LBM::DIR::ZZZ] = mfbbb;
+
+
+ (dist.f[dirE ])[k] = distribution_out.f[VF::LBM::DIR::PZZ];
+ (dist.f[dirW ])[kw] = distribution_out.f[VF::LBM::DIR::MZZ];
+ (dist.f[dirN ])[k] = distribution_out.f[VF::LBM::DIR::ZPZ];
+ (dist.f[dirS ])[ks] = distribution_out.f[VF::LBM::DIR::ZMZ];
+ (dist.f[dirT ])[k] = distribution_out.f[VF::LBM::DIR::ZZP];
+ (dist.f[dirB ])[kb] = distribution_out.f[VF::LBM::DIR::ZZM];
+ (dist.f[dirNE ])[k] = distribution_out.f[VF::LBM::DIR::PPZ];
+ (dist.f[dirSW ])[ksw] = distribution_out.f[VF::LBM::DIR::MMZ];
+ (dist.f[dirSE ])[ks] = distribution_out.f[VF::LBM::DIR::PMZ];
+ (dist.f[dirNW ])[kw] = distribution_out.f[VF::LBM::DIR::MPZ];
+ (dist.f[dirTE ])[k] = distribution_out.f[VF::LBM::DIR::PZP];
+ (dist.f[dirBW ])[kbw] = distribution_out.f[VF::LBM::DIR::MZM];
+ (dist.f[dirBE ])[kb] = distribution_out.f[VF::LBM::DIR::PZM];
+ (dist.f[dirTW ])[kw] = distribution_out.f[VF::LBM::DIR::MZP];
+ (dist.f[dirTN ])[k] = distribution_out.f[VF::LBM::DIR::ZPP];
+ (dist.f[dirBS ])[kbs] = distribution_out.f[VF::LBM::DIR::ZMM];
+ (dist.f[dirBN ])[kb] = distribution_out.f[VF::LBM::DIR::ZPM];
+ (dist.f[dirTS ])[ks] = distribution_out.f[VF::LBM::DIR::ZMP];
+ (dist.f[dirTNE ])[k] = distribution_out.f[VF::LBM::DIR::PPP];
+ (dist.f[dirTNW ])[kw] = distribution_out.f[VF::LBM::DIR::MPP];
+ (dist.f[dirTSE ])[ks] = distribution_out.f[VF::LBM::DIR::PMP];
+ (dist.f[dirTSW ])[ksw] = distribution_out.f[VF::LBM::DIR::MMP];
+ (dist.f[dirBNE ])[kb] = distribution_out.f[VF::LBM::DIR::PPM];
+ (dist.f[dirBNW ])[kbw] = distribution_out.f[VF::LBM::DIR::MPM];
+ (dist.f[dirBSE ])[kbs] = distribution_out.f[VF::LBM::DIR::PMM];
+ (dist.f[dirBSW ])[kbsw] = distribution_out.f[VF::LBM::DIR::MMM];
+ (dist.f[dirREST])[k] = distribution_out.f[VF::LBM::DIR::ZZZ];
+
}
}
////////////////////////////////////////////////////////////////////////////////
-
-////////////////////////////////////////////////////////////////////////////////
-extern "C" __global__ void Cumulant_K17_LBM_Device_Kernel_old(
+extern "C" __global__ void Cumulant_K17_LBM_Device_Kernel_old_origin(
real omega,
uint* typeOfGridNode,
uint* neighborX,
diff --git a/src/lbm/CumulantChimera.h b/src/lbm/CumulantChimera.h
index 90c3d53cd4915e0f5a43b392ac0158387b4c2e8b..690ff227de63ec1ef8872f3f46ee1e9fc117da25 100644
--- a/src/lbm/CumulantChimera.h
+++ b/src/lbm/CumulantChimera.h
@@ -23,16 +23,28 @@ namespace LBM
struct Distribution27
{
- real f[27];
+ real f[27];
- inline __host__ __device__ real getDensity_() const
- {
- return getDensity(f);
- }
+ inline __host__ __device__ real getDensity_() const
+ {
+ return getDensity(f);
+ }
};
-inline __host__ __device__ Distribution27 cumulantChimera(const Distribution27& distribution, real omega, const real* forces)
+
+//////////////////////////////////////////////////////////////////////////
+//! Cumulant K17 Kernel is based on \ref
+//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+//! and \ref
+//! <a href="https://doi.org/10.1016/j.jcp.2017.07.004"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.07.004 ]</b></a>
+//////////////////////////////////////////////////////////////////////////
+inline __host__ __device__ void cumulantChimera(Distribution27& distribution, real omega, real* forces)
{
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Read distributions: style of reading and writing the distributions from/to
+ //! stored arrays dependent on timestep is based on the esoteric twist algorithm
+ //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
+ //!
real mfcbb = distribution.f[DIR::PZZ];
real mfabb = distribution.f[DIR::MZZ];
real mfbcb = distribution.f[DIR::ZPZ];
@@ -66,23 +78,23 @@ inline __host__ __device__ Distribution27 cumulantChimera(const Distribution27&
//! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
//!
real drho =
- ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) +
- (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
- ((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb;
+ ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) +
+ (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
+ ((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb;
real rho = c1o1 + drho;
real OOrho = c1o1 / rho;
real vvx =
- ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
- (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
- (mfcbb - mfabb)) * OOrho;
+ ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+ (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
+ (mfcbb - mfabb)) * OOrho;
real vvy =
- ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
- (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
- (mfbcb - mfbab)) * OOrho;
+ ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+ (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
+ (mfbcb - mfbab)) * OOrho;
real vvz =
- ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
- (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
- (mfbbc - mfbba)) * OOrho;
+ ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+ (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
+ (mfbbc - mfbba)) * OOrho;
////////////////////////////////////////////////////////////////////////////////////
//! - Add half of the acceleration (body force) to the velocity as in Eq. (42) \ref
//! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
@@ -198,17 +210,17 @@ inline __host__ __device__ Distribution27 cumulantChimera(const Distribution27&
////////////////////////////////////////////////////////////
//6.
real CUMccc = mfccc + ((-c4o1 * mfbbb * mfbbb
- - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
- - c4o1 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
- - c2o1 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
- + (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
- + c2o1 * (mfcaa * mfaca * mfaac)
- + c16o1 * mfbba * mfbab * mfabb) * OOrho * OOrho
- - c1o3 * (mfacc + mfcac + mfcca) * OOrho
- - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
- + (c2o1 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
- + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) * OOrho * OOrho * c2o3
- + c1o27*((drho * drho - drho) * OOrho * OOrho));
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - c4o1 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+ - c2o1 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
+ + (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + c2o1 * (mfcaa * mfaca * mfaac)
+ + c16o1 * mfbba * mfbab * mfabb) * OOrho * OOrho
+ - c1o3 * (mfacc + mfcac + mfcca) * OOrho
+ - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
+ + (c2o1 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
+ + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) * OOrho * OOrho * c2o3
+ + c1o27*((drho * drho - drho) * OOrho * OOrho));
////////////////////////////////////////////////////////////////////////////////////
//! - Compute linear combinations of second and third order cumulants
//!
@@ -339,17 +351,17 @@ inline __host__ __device__ Distribution27 cumulantChimera(const Distribution27&
//////////////////////////////////////////////////////////////////////////
//6.
mfccc = CUMccc - ((-c4o1 * mfbbb * mfbbb
- - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
- - c4o1 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
- - c2o1 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
- + (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
- + c2o1 * (mfcaa * mfaca * mfaac)
- + c16o1 * mfbba * mfbab * mfabb) * OOrho * OOrho
- - c1o3 * (mfacc + mfcac + mfcca) * OOrho
- - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
- + (c2o1 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
- + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) * OOrho * OOrho * c2o3
- + c1o27*((drho * drho - drho) * OOrho * OOrho));
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - c4o1 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+ - c2o1 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
+ + (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + c2o1 * (mfcaa * mfaca * mfaac)
+ + c16o1 * mfbba * mfbab * mfabb) * OOrho * OOrho
+ - c1o3 * (mfacc + mfcac + mfcca) * OOrho
+ - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
+ + (c2o1 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
+ + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) * OOrho * OOrho * c2o3
+ + c1o27*((drho * drho - drho) * OOrho * OOrho));
////////////////////////////////////////////////////////////////////////////////////
//! - Add acceleration (body force) to first order cumulants according to Eq. (85)-(87) in
//! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
@@ -397,39 +409,42 @@ inline __host__ __device__ Distribution27 cumulantChimera(const Distribution27&
VF::LBM::backwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9o1, c1o9);
VF::LBM::backwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36o1, c1o36);
- Distribution27 distribution_calculated;
-
- distribution_calculated.f[DIR::PZZ] = mfcbb;
- distribution_calculated.f[DIR::MZZ] = mfabb;
- distribution_calculated.f[DIR::ZPZ] = mfbcb;
- distribution_calculated.f[DIR::ZMZ] = mfbab;
- distribution_calculated.f[DIR::ZZP] = mfbbc;
- distribution_calculated.f[DIR::ZZM] = mfbba;
- distribution_calculated.f[DIR::PPZ] = mfccb;
- distribution_calculated.f[DIR::MMZ] = mfaab;
- distribution_calculated.f[DIR::PMZ] = mfcab;
- distribution_calculated.f[DIR::MPZ] = mfacb;
- distribution_calculated.f[DIR::PZP] = mfcbc;
- distribution_calculated.f[DIR::MZM] = mfaba;
- distribution_calculated.f[DIR::PZM] = mfcba;
- distribution_calculated.f[DIR::MZP] = mfabc;
- distribution_calculated.f[DIR::ZPP] = mfbcc;
- distribution_calculated.f[DIR::ZMM] = mfbaa;
- distribution_calculated.f[DIR::ZPM] = mfbca;
- distribution_calculated.f[DIR::ZMP] = mfbac;
- distribution_calculated.f[DIR::PPP] = mfccc;
- distribution_calculated.f[DIR::MPP] = mfacc;
- distribution_calculated.f[DIR::PMP] = mfcac;
- distribution_calculated.f[DIR::MMP] = mfaac;
- distribution_calculated.f[DIR::PPM] = mfcca;
- distribution_calculated.f[DIR::MPM] = mfaca;
- distribution_calculated.f[DIR::PMM] = mfcaa;
- distribution_calculated.f[DIR::MMM] = mfaaa;
- distribution_calculated.f[DIR::ZZZ] = mfbbb;
- return distribution_calculated;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Write distributions: style of reading and writing the distributions from/to
+ //! stored arrays dependent on timestep is based on the esoteric twist algorithm
+ //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
+ //!
+ distribution.f[VF::LBM::DIR::MZZ] = mfcbb;
+ distribution.f[VF::LBM::DIR::PZZ] = mfabb;
+ distribution.f[VF::LBM::DIR::ZMZ] = mfbcb;
+ distribution.f[VF::LBM::DIR::ZPZ] = mfbab;
+ distribution.f[VF::LBM::DIR::ZZM] = mfbbc;
+ distribution.f[VF::LBM::DIR::ZZP] = mfbba;
+ distribution.f[VF::LBM::DIR::MMZ] = mfccb;
+ distribution.f[VF::LBM::DIR::PPZ] = mfaab;
+ distribution.f[VF::LBM::DIR::MPZ] = mfcab;
+ distribution.f[VF::LBM::DIR::PMZ] = mfacb;
+ distribution.f[VF::LBM::DIR::MZM] = mfcbc;
+ distribution.f[VF::LBM::DIR::PZP] = mfaba;
+ distribution.f[VF::LBM::DIR::MZP] = mfcba;
+ distribution.f[VF::LBM::DIR::PZM] = mfabc;
+ distribution.f[VF::LBM::DIR::ZMM] = mfbcc;
+ distribution.f[VF::LBM::DIR::ZPP] = mfbaa;
+ distribution.f[VF::LBM::DIR::ZMP] = mfbca;
+ distribution.f[VF::LBM::DIR::ZPM] = mfbac;
+ distribution.f[VF::LBM::DIR::MMM] = mfccc;
+ distribution.f[VF::LBM::DIR::PMM] = mfacc;
+ distribution.f[VF::LBM::DIR::MPM] = mfcac;
+ distribution.f[VF::LBM::DIR::PPM] = mfaac;
+ distribution.f[VF::LBM::DIR::MMP] = mfcca;
+ distribution.f[VF::LBM::DIR::PMP] = mfaca;
+ distribution.f[VF::LBM::DIR::MPP] = mfcaa;
+ distribution.f[VF::LBM::DIR::PPP] = mfaaa;
+ distribution.f[VF::LBM::DIR::ZZZ] = mfbbb;
}
+
}
}
#endif