diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim.cu b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim.cu
index fb39d57db2101e1edecf1a7707b5bc9c306c6b1c..0980050dcb43f84b5a7da42dd2c9e4cef1796833 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim.cu
@@ -15,12 +15,6 @@ void TurbulentViscosityCumulantK17CompChim::run()
TurbulenceModel turbulenceModel = para->getTurbulenceModel();
switch(para->getTurbulenceModel())
{
- case TurbulenceModel::None:
- LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::None > <<< grid.grid, grid.threads >>>( para->getParD(level)->omega, para->getParD(level)->typeOfGridNode, para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ, para->getParD(level)->distributions.f[0],
- para->getParD(level)->rho, para->getParD(level)->velocityX, para->getParD(level)->velocityY, para->getParD(level)->velocityZ, para->getParD(level)->turbViscosity,para->getSGSConstant(),
- (unsigned long)para->getParD(level)->numberOfNodes, level, para->getIsBodyForce(), para->getForcesDev(), para->getParD(level)->forceX_SP, para->getParD(level)->forceY_SP,
- para->getParD(level)->forceZ_SP,para->getQuadricLimitersDev(), para->getParD(level)->isEvenTimestep);
- break;
case TurbulenceModel::AMD:
LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::AMD > <<< grid.grid, grid.threads >>>( para->getParD(level)->omega, para->getParD(level)->typeOfGridNode, para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ, para->getParD(level)->distributions.f[0],
para->getParD(level)->rho, para->getParD(level)->velocityX, para->getParD(level)->velocityY, para->getParD(level)->velocityZ, para->getParD(level)->turbViscosity,para->getSGSConstant(),
@@ -39,8 +33,11 @@ void TurbulentViscosityCumulantK17CompChim::run()
(unsigned long)para->getParD(level)->numberOfNodes, level, para->getIsBodyForce(), para->getForcesDev(), para->getParD(level)->forceX_SP, para->getParD(level)->forceY_SP,
para->getParD(level)->forceZ_SP,para->getQuadricLimitersDev(), para->getParD(level)->isEvenTimestep);
break;
+ case TurbulenceModel::None:
+ throw std::runtime_error("TurbulentViscosityCumulantK17CompChim cannot be use without turbulence Model (TurbulenceModel::None)!");
+ break;
default:
- throw std::runtime_error("TurbulentViscosityCumulantK17CompChim: Invalid turbulence mpdel ");
+ throw std::runtime_error("TurbulentViscosityCumulantK17CompChim: Invalid turbulence model!");
break;
}
getLastCudaError("LB_Kernel_TurbulentViscosityCumulantK17CompChim execution failed");
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cu b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cu
index af7a67ce3f0b0cd9c4c6279d0728ab838c76fd4a..afab2fba994738fdc7fb509a7611df52be722abf 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/TurbulentViscosityKernels/FluidFlow/Compressible/CumulantK17chim/TurbulentViscosityCumulantK17CompChim_Device.cu
@@ -42,6 +42,8 @@
#include "LBM/LB.h"
#include "lbm/constants/D3Q27.h"
#include <lbm/constants/NumericConstants.h>
+#include "Kernel/Utilities/DistributionHelper.cuh"
+
#include "GPU/TurbulentViscosityInlines.cuh"
using namespace vf::lbm::constant;
@@ -85,144 +87,112 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
////////////////////////////////////////////////////////////////////////////////
//! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
//!
- 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;
-
- const unsigned k = nx * (ny * z + y) + x;
+ const unsigned k_000 = vf::gpu::getNodeIndex();
//////////////////////////////////////////////////////////////////////////
// run for all indices in size_Mat and fluid nodes
- if ((k < size_Mat) && (typeOfGridNode[k] == GEO_FLUID)) {
+ if ((k_000 < size_Mat) && (typeOfGridNode[k_000] == GEO_FLUID)) {
//////////////////////////////////////////////////////////////////////////
//! - Read distributions: style of reading and writing the distributions from/to stored arrays dependent on
//! timestep is based on the esoteric twist algorithm \ref <a
//! href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017),
//! DOI:10.3390/computation5020019 ]</b></a>
//!
- Distributions27 dist;
- if (isEvenTimestep) {
- dist.f[DIR_P00] = &distributions[DIR_P00 * size_Mat];
- dist.f[DIR_M00] = &distributions[DIR_M00 * size_Mat];
- dist.f[DIR_0P0] = &distributions[DIR_0P0 * size_Mat];
- dist.f[DIR_0M0] = &distributions[DIR_0M0 * size_Mat];
- dist.f[DIR_00P] = &distributions[DIR_00P * size_Mat];
- dist.f[DIR_00M] = &distributions[DIR_00M * size_Mat];
- dist.f[DIR_PP0] = &distributions[DIR_PP0 * size_Mat];
- dist.f[DIR_MM0] = &distributions[DIR_MM0 * size_Mat];
- dist.f[DIR_PM0] = &distributions[DIR_PM0 * size_Mat];
- dist.f[DIR_MP0] = &distributions[DIR_MP0 * size_Mat];
- dist.f[DIR_P0P] = &distributions[DIR_P0P * size_Mat];
- dist.f[DIR_M0M] = &distributions[DIR_M0M * size_Mat];
- dist.f[DIR_P0M] = &distributions[DIR_P0M * size_Mat];
- dist.f[DIR_M0P] = &distributions[DIR_M0P * size_Mat];
- dist.f[DIR_0PP] = &distributions[DIR_0PP * size_Mat];
- dist.f[DIR_0MM] = &distributions[DIR_0MM * size_Mat];
- dist.f[DIR_0PM] = &distributions[DIR_0PM * size_Mat];
- dist.f[DIR_0MP] = &distributions[DIR_0MP * size_Mat];
- dist.f[DIR_000] = &distributions[DIR_000 * size_Mat];
- dist.f[DIR_PPP] = &distributions[DIR_PPP * size_Mat];
- dist.f[DIR_MMP] = &distributions[DIR_MMP * size_Mat];
- dist.f[DIR_PMP] = &distributions[DIR_PMP * size_Mat];
- dist.f[DIR_MPP] = &distributions[DIR_MPP * size_Mat];
- dist.f[DIR_PPM] = &distributions[DIR_PPM * size_Mat];
- dist.f[DIR_MMM] = &distributions[DIR_MMM * size_Mat];
- dist.f[DIR_PMM] = &distributions[DIR_PMM * size_Mat];
- dist.f[DIR_MPM] = &distributions[DIR_MPM * size_Mat];
- } else {
- dist.f[DIR_M00] = &distributions[DIR_P00 * size_Mat];
- dist.f[DIR_P00] = &distributions[DIR_M00 * size_Mat];
- dist.f[DIR_0M0] = &distributions[DIR_0P0 * size_Mat];
- dist.f[DIR_0P0] = &distributions[DIR_0M0 * size_Mat];
- dist.f[DIR_00M] = &distributions[DIR_00P * size_Mat];
- dist.f[DIR_00P] = &distributions[DIR_00M * size_Mat];
- dist.f[DIR_MM0] = &distributions[DIR_PP0 * size_Mat];
- dist.f[DIR_PP0] = &distributions[DIR_MM0 * size_Mat];
- dist.f[DIR_MP0] = &distributions[DIR_PM0 * size_Mat];
- dist.f[DIR_PM0] = &distributions[DIR_MP0 * size_Mat];
- dist.f[DIR_M0M] = &distributions[DIR_P0P * size_Mat];
- dist.f[DIR_P0P] = &distributions[DIR_M0M * size_Mat];
- dist.f[DIR_M0P] = &distributions[DIR_P0M * size_Mat];
- dist.f[DIR_P0M] = &distributions[DIR_M0P * size_Mat];
- dist.f[DIR_0MM] = &distributions[DIR_0PP * size_Mat];
- dist.f[DIR_0PP] = &distributions[DIR_0MM * size_Mat];
- dist.f[DIR_0MP] = &distributions[DIR_0PM * size_Mat];
- dist.f[DIR_0PM] = &distributions[DIR_0MP * size_Mat];
- dist.f[DIR_000] = &distributions[DIR_000 * size_Mat];
- dist.f[DIR_MMM] = &distributions[DIR_PPP * size_Mat];
- dist.f[DIR_PPM] = &distributions[DIR_MMP * size_Mat];
- dist.f[DIR_MPM] = &distributions[DIR_PMP * size_Mat];
- dist.f[DIR_PMM] = &distributions[DIR_MPP * size_Mat];
- dist.f[DIR_MMP] = &distributions[DIR_PPM * size_Mat];
- dist.f[DIR_PPP] = &distributions[DIR_MMM * size_Mat];
- dist.f[DIR_MPP] = &distributions[DIR_PMM * size_Mat];
- dist.f[DIR_PMP] = &distributions[DIR_MPM * size_Mat];
- }
+ Distributions27 dist = vf::gpu::getDistributionReferences27(distributions, size_Mat, isEvenTimestep);
+
////////////////////////////////////////////////////////////////////////////////
//! - 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];
+ uint k_M00 = neighborX[k_000];
+ uint k_0M0 = neighborY[k_000];
+ uint k_00M = neighborZ[k_000];
+ uint k_MM0 = neighborY[k_M00];
+ uint k_M0M = neighborZ[k_M00];
+ uint k_0MM = neighborZ[k_0M0];
+ uint k_MMM = neighborZ[k_MM0];
////////////////////////////////////////////////////////////////////////////////////
//! - Set local distributions
//!
- real mfcbb = (dist.f[DIR_P00])[k];
- real mfabb = (dist.f[DIR_M00])[kw];
- real mfbcb = (dist.f[DIR_0P0])[k];
- real mfbab = (dist.f[DIR_0M0])[ks];
- real mfbbc = (dist.f[DIR_00P])[k];
- real mfbba = (dist.f[DIR_00M])[kb];
- real mfccb = (dist.f[DIR_PP0])[k];
- real mfaab = (dist.f[DIR_MM0])[ksw];
- real mfcab = (dist.f[DIR_PM0])[ks];
- real mfacb = (dist.f[DIR_MP0])[kw];
- real mfcbc = (dist.f[DIR_P0P])[k];
- real mfaba = (dist.f[DIR_M0M])[kbw];
- real mfcba = (dist.f[DIR_P0M])[kb];
- real mfabc = (dist.f[DIR_M0P])[kw];
- real mfbcc = (dist.f[DIR_0PP])[k];
- real mfbaa = (dist.f[DIR_0MM])[kbs];
- real mfbca = (dist.f[DIR_0PM])[kb];
- real mfbac = (dist.f[DIR_0MP])[ks];
- real mfbbb = (dist.f[DIR_000])[k];
- real mfccc = (dist.f[DIR_PPP])[k];
- real mfaac = (dist.f[DIR_MMP])[ksw];
- real mfcac = (dist.f[DIR_PMP])[ks];
- real mfacc = (dist.f[DIR_MPP])[kw];
- real mfcca = (dist.f[DIR_PPM])[kb];
- real mfaaa = (dist.f[DIR_MMM])[kbsw];
- real mfcaa = (dist.f[DIR_PMM])[kbs];
- real mfaca = (dist.f[DIR_MPM])[kbw];
+ real f_000 = (dist.f[DIR_000])[k_000];
+ real f_P00 = (dist.f[DIR_P00])[k_000];
+ real f_M00 = (dist.f[DIR_M00])[k_M00];
+ real f_0P0 = (dist.f[DIR_0P0])[k_000];
+ real f_0M0 = (dist.f[DIR_0M0])[k_0M0];
+ real f_00P = (dist.f[DIR_00P])[k_000];
+ real f_00M = (dist.f[DIR_00M])[k_00M];
+ real f_PP0 = (dist.f[DIR_PP0])[k_000];
+ real f_MM0 = (dist.f[DIR_MM0])[k_MM0];
+ real f_PM0 = (dist.f[DIR_PM0])[k_0M0];
+ real f_MP0 = (dist.f[DIR_MP0])[k_M00];
+ real f_P0P = (dist.f[DIR_P0P])[k_000];
+ real f_M0M = (dist.f[DIR_M0M])[k_M0M];
+ real f_P0M = (dist.f[DIR_P0M])[k_00M];
+ real f_M0P = (dist.f[DIR_M0P])[k_M00];
+ real f_0PP = (dist.f[DIR_0PP])[k_000];
+ real f_0MM = (dist.f[DIR_0MM])[k_0MM];
+ real f_0PM = (dist.f[DIR_0PM])[k_00M];
+ real f_0MP = (dist.f[DIR_0MP])[k_0M0];
+ real f_PPP = (dist.f[DIR_PPP])[k_000];
+ real f_MPP = (dist.f[DIR_MPP])[k_M00];
+ real f_PMP = (dist.f[DIR_PMP])[k_0M0];
+ real f_MMP = (dist.f[DIR_MMP])[k_MM0];
+ real f_PPM = (dist.f[DIR_PPM])[k_00M];
+ real f_MPM = (dist.f[DIR_MPM])[k_M0M];
+ real f_PMM = (dist.f[DIR_PMM])[k_0MM];
+ real f_MMM = (dist.f[DIR_MMM])[k_MMM];
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Define aliases to use the same variable for the moments (m's):
+ //!
+ real& m_111 = f_000;
+ real& m_211 = f_P00;
+ real& m_011 = f_M00;
+ real& m_121 = f_0P0;
+ real& m_101 = f_0M0;
+ real& m_112 = f_00P;
+ real& m_110 = f_00M;
+ real& m_221 = f_PP0;
+ real& m_001 = f_MM0;
+ real& m_201 = f_PM0;
+ real& m_021 = f_MP0;
+ real& m_212 = f_P0P;
+ real& m_010 = f_M0M;
+ real& m_210 = f_P0M;
+ real& m_012 = f_M0P;
+ real& m_122 = f_0PP;
+ real& m_100 = f_0MM;
+ real& m_120 = f_0PM;
+ real& m_102 = f_0MP;
+ real& m_222 = f_PPP;
+ real& m_022 = f_MPP;
+ real& m_202 = f_PMP;
+ real& m_002 = f_MMP;
+ real& m_220 = f_PPM;
+ real& m_020 = f_MPM;
+ real& m_200 = f_PMM;
+ real& m_000 = f_MMM;
+
//////////////////////////////////////////////////////(unsigned long)//////////////////////////////
//! - 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 rrho = c1o1 + drho;
- real OOrho = c1o1 / rrho;
-
- 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;
+ real drho = ((((f_PPP + f_MMM) + (f_MPM + f_PMP)) + ((f_MPP + f_PMM) + (f_MMP + f_PPM))) +
+ (((f_0MP + f_0PM) + (f_0MM + f_0PP)) + ((f_M0P + f_P0M) + (f_M0M + f_P0P)) +
+ ((f_MP0 + f_PM0) + (f_MM0 + f_PP0))) +
+ ((f_M00 + f_P00) + (f_0M0 + f_0P0) + (f_00M + f_00P))) +
+ f_000;
+
+ real oneOverRho = c1o1 / (c1o1 + drho);
+
+ real vvx = ((((f_PPP - f_MMM) + (f_PMP - f_MPM)) + ((f_PMM - f_MPP) + (f_PPM - f_MMP))) +
+ (((f_P0M - f_M0P) + (f_P0P - f_M0M)) + ((f_PM0 - f_MP0) + (f_PP0 - f_MM0))) + (f_P00 - f_M00)) *
+ oneOverRho;
+ real vvy = ((((f_PPP - f_MMM) + (f_MPM - f_PMP)) + ((f_MPP - f_PMM) + (f_PPM - f_MMP))) +
+ (((f_0PM - f_0MP) + (f_0PP - f_0MM)) + ((f_MP0 - f_PM0) + (f_PP0 - f_MM0))) + (f_0P0 - f_0M0)) *
+ oneOverRho;
+ real vvz = ((((f_PPP - f_MMM) + (f_PMP - f_MPM)) + ((f_MPP - f_PMM) + (f_MMP - f_PPM))) +
+ (((f_0MP - f_0PM) + (f_0PP - f_0MM)) + ((f_M0P - f_P0M) + (f_P0P - f_M0M))) + (f_00P - f_00M)) *
+ oneOverRho;
////////////////////////////////////////////////////////////////////////////////////
//! - Add half of the acceleration (body force) to the velocity as in Eq. (42) \ref
@@ -239,9 +209,9 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
real fz = forces[2];
if( bodyForce ){
- fx += bodyForceX[k];
- fy += bodyForceY[k];
- fz += bodyForceZ[k];
+ fx += bodyForceX[k_000];
+ fy += bodyForceY[k_000];
+ fz += bodyForceZ[k_000];
real vx = vvx;
real vy = vvy;
@@ -268,9 +238,9 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//!> differ by several orders of magnitude.
//!> \note 16/05/2022: Testing, still ongoing!
//!
- bodyForceX[k] = (acc_x-(vvx-vx))*factor*c2o1;
- bodyForceY[k] = (acc_y-(vvy-vy))*factor*c2o1;
- bodyForceZ[k] = (acc_z-(vvz-vz))*factor*c2o1;
+ bodyForceX[k_000] = (acc_x-(vvx-vx))*factor*c2o1;
+ bodyForceY[k_000] = (acc_y-(vvy-vy))*factor*c2o1;
+ bodyForceZ[k_000] = (acc_z-(vvz-vz))*factor*c2o1;
}
else{
vvx += fx * c1o2 / factor;
@@ -289,10 +259,9 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//! 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 = quadricLimiters[0];
- real qudricLimitM = quadricLimiters[1];
- real qudricLimitD = quadricLimiters[2];
+ real quadricLimitP = quadricLimiters[0];
+ real quadricLimitM = quadricLimiters[1];
+ real quadricLimitD = quadricLimiters[2];
////////////////////////////////////////////////////////////////////////////////////
//! - 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),
@@ -302,39 +271,39 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//!
////////////////////////////////////////////////////////////////////////////////////
// Z - Dir
- forwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36o1, c1o36);
- forwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9o1, c1o9);
- forwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36o1, c1o36);
- forwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9o1, c1o9);
- forwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
- forwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9o1, c1o9);
- forwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36o1, c1o36);
- forwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9o1, c1o9);
- forwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36o1, c1o36);
+ forwardInverseChimeraWithK(f_MMM, f_MM0, f_MMP, vvz, vz2, c36o1, c1o36);
+ forwardInverseChimeraWithK(f_M0M, f_M00, f_M0P, vvz, vz2, c9o1, c1o9);
+ forwardInverseChimeraWithK(f_MPM, f_MP0, f_MPP, vvz, vz2, c36o1, c1o36);
+ forwardInverseChimeraWithK(f_0MM, f_0M0, f_0MP, vvz, vz2, c9o1, c1o9);
+ forwardInverseChimeraWithK(f_00M, f_000, f_00P, vvz, vz2, c9o4, c4o9);
+ forwardInverseChimeraWithK(f_0PM, f_0P0, f_0PP, vvz, vz2, c9o1, c1o9);
+ forwardInverseChimeraWithK(f_PMM, f_PM0, f_PMP, vvz, vz2, c36o1, c1o36);
+ forwardInverseChimeraWithK(f_P0M, f_P00, f_P0P, vvz, vz2, c9o1, c1o9);
+ forwardInverseChimeraWithK(f_PPM, f_PP0, f_PPP, vvz, vz2, c36o1, c1o36);
////////////////////////////////////////////////////////////////////////////////////
// Y - Dir
- forwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6o1, c1o6);
- forwardChimera(mfaab, mfabb, mfacb, vvy, vy2);
- forwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18o1, c1o18);
- forwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
- forwardChimera(mfbab, mfbbb, mfbcb, vvy, vy2);
- forwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
- forwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6o1, c1o6);
- forwardChimera(mfcab, mfcbb, mfccb, vvy, vy2);
- forwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18o1, c1o18);
+ forwardInverseChimeraWithK(f_MMM, f_M0M, f_MPM, vvy, vy2, c6o1, c1o6);
+ forwardChimera( f_MM0, f_M00, f_MP0, vvy, vy2);
+ forwardInverseChimeraWithK(f_MMP, f_M0P, f_MPP, vvy, vy2, c18o1, c1o18);
+ forwardInverseChimeraWithK(f_0MM, f_00M, f_0PM, vvy, vy2, c3o2, c2o3);
+ forwardChimera( f_0M0, f_000, f_0P0, vvy, vy2);
+ forwardInverseChimeraWithK(f_0MP, f_00P, f_0PP, vvy, vy2, c9o2, c2o9);
+ forwardInverseChimeraWithK(f_PMM, f_P0M, f_PPM, vvy, vy2, c6o1, c1o6);
+ forwardChimera( f_PM0, f_P00, f_PP0, vvy, vy2);
+ forwardInverseChimeraWithK(f_PMP, f_P0P, f_PPP, vvy, vy2, c18o1, c1o18);
////////////////////////////////////////////////////////////////////////////////////
// X - Dir
- forwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1o1, c1o1);
- forwardChimera(mfaba, mfbba, mfcba, vvx, vx2);
- forwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3o1, c1o3);
- forwardChimera(mfaab, mfbab, mfcab, vvx, vx2);
- forwardChimera(mfabb, mfbbb, mfcbb, vvx, vx2);
- forwardChimera(mfacb, mfbcb, mfccb, vvx, vx2);
- forwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3o1, c1o3);
- forwardChimera(mfabc, mfbbc, mfcbc, vvx, vx2);
- forwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c3o1, c1o9);
+ forwardInverseChimeraWithK(f_MMM, f_0MM, f_PMM, vvx, vx2, c1o1, c1o1);
+ forwardChimera( f_M0M, f_00M, f_P0M, vvx, vx2);
+ forwardInverseChimeraWithK(f_MPM, f_0PM, f_PPM, vvx, vx2, c3o1, c1o3);
+ forwardChimera( f_MM0, f_0M0, f_PM0, vvx, vx2);
+ forwardChimera( f_M00, f_000, f_P00, vvx, vx2);
+ forwardChimera( f_MP0, f_0P0, f_PP0, vvx, vx2);
+ forwardInverseChimeraWithK(f_MMP, f_0MP, f_PMP, vvx, vx2, c3o1, c1o3);
+ forwardChimera( f_M0P, f_00P, f_P0P, vvx, vx2);
+ forwardInverseChimeraWithK(f_MPP, f_0PP, f_PPP, vvx, vx2, c3o1, c1o9);
////////////////////////////////////////////////////////////////////////////////////
//! - Setting relaxation rates for non-hydrodynamic cumulants (default values). Variable names and equations
@@ -357,7 +326,7 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
////////////////////////////////////////////////////////////////////////////////////
//! - Calculate modified omega with turbulent viscosity
//!
- real omega = omega_in / (c1o1 + c3o1*omega_in*turbulentViscosity[k]);
+ real omega = omega_in / (c1o1 + c3o1*omega_in*turbulentViscosity[k_000]);
////////////////////////////////////////////////////////////
// 2.
real OxxPyyPzz = c1o1;
@@ -394,63 +363,60 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//!
////////////////////////////////////////////////////////////
// 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));
+ real c_211 = m_211 - ((m_200 + c1o3) * m_011 + c2o1 * m_110 * m_101) * oneOverRho;
+ real c_121 = m_121 - ((m_020 + c1o3) * m_101 + c2o1 * m_110 * m_011) * oneOverRho;
+ real c_112 = m_112 - ((m_002 + c1o3) * m_110 + c2o1 * m_101 * m_011) * oneOverRho;
+
+ real c_220 = m_220 - (((m_200 * m_020 + c2o1 * m_110 * m_110) + c1o3 * (m_200 + m_020)) * oneOverRho - c1o9 * (drho * oneOverRho));
+ real c_202 = m_202 - (((m_200 * m_002 + c2o1 * m_101 * m_101) + c1o3 * (m_200 + m_002)) * oneOverRho - c1o9 * (drho * oneOverRho));
+ real c_022 = m_022 - (((m_002 * m_020 + c2o1 * m_011 * m_011) + c1o3 * (m_002 + m_020)) * oneOverRho - c1o9 * (drho * oneOverRho));
////////////////////////////////////////////////////////////
// 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;
+ real c_122 =
+ m_122 - ((m_002 * m_120 + m_020 * m_102 + c4o1 * m_011 * m_111 + c2o1 * (m_101 * m_021 + m_110 * m_012)) +
+ c1o3 * (m_120 + m_102)) *
+ oneOverRho;
+ real c_212 =
+ m_212 - ((m_002 * m_210 + m_200 * m_012 + c4o1 * m_101 * m_111 + c2o1 * (m_011 * m_201 + m_110 * m_102)) +
+ c1o3 * (m_210 + m_012)) *
+ oneOverRho;
+ real c_221 =
+ m_221 - ((m_200 * m_021 + m_020 * m_201 + c4o1 * m_110 * m_111 + c2o1 * (m_101 * m_120 + m_011 * m_210)) +
+ c1o3 * (m_021 + m_201)) *
+ oneOverRho;
////////////////////////////////////////////////////////////
// 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));
+ real c_222 = m_222 + ((-c4o1 * m_111 * m_111 - (m_200 * m_022 + m_020 * m_202 + m_002 * m_220) -
+ c4o1 * (m_011 * m_211 + m_101 * m_121 + m_110 * m_112) -
+ c2o1 * (m_120 * m_102 + m_210 * m_012 + m_201 * m_021)) *
+ oneOverRho +
+ (c4o1 * (m_101 * m_101 * m_020 + m_011 * m_011 * m_200 + m_110 * m_110 * m_002) +
+ c2o1 * (m_200 * m_020 * m_002) + c16o1 * m_110 * m_101 * m_011) *
+ oneOverRho * oneOverRho -
+ c1o3 * (m_022 + m_202 + m_220) * oneOverRho - c1o9 * (m_200 + m_020 + m_002) * oneOverRho +
+ (c2o1 * (m_101 * m_101 + m_011 * m_011 + m_110 * m_110) +
+ (m_002 * m_020 + m_002 * m_200 + m_020 * m_200) + c1o3 * (m_002 + m_020 + m_200)) *
+ oneOverRho * oneOverRho * c2o3 +
+ c1o27 * ((drho * drho - drho) * oneOverRho * oneOverRho));
////////////////////////////////////////////////////////////////////////////////////
//! - Compute linear combinations of second and third order cumulants
//!
////////////////////////////////////////////////////////////
// 2.
- real mxxPyyPzz = mfcaa + mfaca + mfaac;
- real mxxMyy = mfcaa - mfaca;
- real mxxMzz = mfcaa - mfaac;
+ real mxxPyyPzz = m_200 + m_020 + m_002;
+ real mxxMyy = m_200 - m_020;
+ real mxxMzz = m_200 - m_002;
////////////////////////////////////////////////////////////
// 3.
- real mxxyPyzz = mfcba + mfabc;
- real mxxyMyzz = mfcba - mfabc;
+ real mxxyPyzz = m_210 + m_012;
+ real mxxyMyzz = m_210 - m_012;
- real mxxzPyyz = mfcab + mfacb;
- real mxxzMyyz = mfcab - mfacb;
+ real mxxzPyyz = m_201 + m_021;
+ real mxxzMyyz = m_201 - m_021;
- real mxyyPxzz = mfbca + mfbac;
- real mxyyMxzz = mfbca - mfbac;
+ real mxyyPxzz = m_120 + m_102;
+ real mxyyMxzz = m_120 - m_102;
////////////////////////////////////////////////////////////////////////////////////
// incl. correction
@@ -462,24 +428,23 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//! 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 Dxy = -c3o1 * omega * m_110;
+ real Dxz = -c3o1 * omega * m_101;
+ real Dyz = -c3o1 * omega * m_011;
+ real dxux = c1o2 * (-omega) * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (m_000 - mxxPyyPzz);
real dyuy = dxux + omega * c3o2 * mxxMyy;
real dzuz = dxux + omega * c3o2 * mxxMzz;
////////////////////////////////////////////////////////////////////////////////////
switch (turbulenceModel)
{
- case TurbulenceModel::None:
case TurbulenceModel::AMD: //AMD is computed in separate kernel
break;
case TurbulenceModel::Smagorinsky:
- turbulentViscosity[k] = calcTurbulentViscositySmagorinsky(SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
+ turbulentViscosity[k_000] = calcTurbulentViscositySmagorinsky(SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
break;
case TurbulenceModel::QR:
- turbulentViscosity[k] = calcTurbulentViscosityQR(SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
+ turbulentViscosity[k_000] = calcTurbulentViscosityQR(SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
break;
default:
break;
@@ -489,8 +454,7 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//! <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);
+ mxxPyyPzz += OxxPyyPzz * (m_000 - 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);
@@ -500,9 +464,9 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
// mxxMyy += -(-omega) * (-mxxMyy);
// mxxMzz += -(-omega) * (-mxxMzz);
//////////////////////////////////////////////////////////////////////////
- mfabb += omega * (-mfabb);
- mfbab += omega * (-mfbab);
- mfbba += omega * (-mfbba);
+ m_011 += omega * (-m_011);
+ m_101 += omega * (-m_101);
+ m_110 += omega * (-m_110);
////////////////////////////////////////////////////////////////////////////////////
// relax
@@ -512,19 +476,19 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//! <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);
+ real wadjust = Oxyz + (c1o1 - Oxyz) * abs(m_111) / (abs(m_111) + quadricLimitD);
+ m_111 += wadjust * (-m_111);
+ wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * abs(mxxyPyzz) / (abs(mxxyPyzz) + quadricLimitP);
mxxyPyzz += wadjust * (-mxxyPyzz);
- wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * abs(mxxyMyzz) / (abs(mxxyMyzz) + qudricLimitM);
+ wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * abs(mxxyMyzz) / (abs(mxxyMyzz) + quadricLimitM);
mxxyMyzz += wadjust * (-mxxyMyzz);
- wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * abs(mxxzPyyz) / (abs(mxxzPyyz) + qudricLimitP);
+ wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * abs(mxxzPyyz) / (abs(mxxzPyyz) + quadricLimitP);
mxxzPyyz += wadjust * (-mxxzPyyz);
- wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * abs(mxxzMyyz) / (abs(mxxzMyyz) + qudricLimitM);
+ wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * abs(mxxzMyyz) / (abs(mxxzMyyz) + quadricLimitM);
mxxzMyyz += wadjust * (-mxxzMyyz);
- wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * abs(mxyyPxzz) / (abs(mxyyPxzz) + qudricLimitP);
+ wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * abs(mxyyPxzz) / (abs(mxyyPxzz) + quadricLimitP);
mxyyPxzz += wadjust * (-mxyyPxzz);
- wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * abs(mxyyMxzz) / (abs(mxyyMxzz) + qudricLimitM);
+ wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * abs(mxyyMxzz) / (abs(mxyyMxzz) + quadricLimitM);
mxyyMxzz += wadjust * (-mxyyMxzz);
//////////////////////////////////////////////////////////////////////////
// no limiter
@@ -539,16 +503,16 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
////////////////////////////////////////////////////////////////////////////////////
//! - 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;
+ m_200 = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz);
+ m_020 = c1o3 * (-c2o1 * mxxMyy + mxxMzz + mxxPyyPzz);
+ m_002 = c1o3 * (mxxMyy - c2o1 * mxxMzz + mxxPyyPzz);
+
+ m_210 = ( mxxyMyzz + mxxyPyzz) * c1o2;
+ m_012 = (-mxxyMyzz + mxxyPyzz) * c1o2;
+ m_201 = ( mxxzMyyz + mxxzPyyz) * c1o2;
+ m_021 = (-mxxzMyyz + mxxzPyyz) * c1o2;
+ m_120 = ( mxyyMxzz + mxyyPxzz) * c1o2;
+ m_102 = (-mxyyMxzz + mxyyPxzz) * c1o2;
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
@@ -558,22 +522,23 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//! 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 * factorA + (c1o1 - O4) * (CUMacc);
- CUMcac = -O4 * (c1o1 / omega - c1o2) * (dxux + dzuz) * c2o3 * factorA + (c1o1 - O4) * (CUMcac);
- CUMcca = -O4 * (c1o1 / omega - c1o2) * (dyuy + dxux) * c2o3 * factorA + (c1o1 - O4) * (CUMcca);
- CUMbbc = -O4 * (c1o1 / omega - c1o2) * Dxy * c1o3 * factorB + (c1o1 - O4) * (CUMbbc);
- CUMbcb = -O4 * (c1o1 / omega - c1o2) * Dxz * c1o3 * factorB + (c1o1 - O4) * (CUMbcb);
- CUMcbb = -O4 * (c1o1 / omega - c1o2) * Dyz * c1o3 * factorB + (c1o1 - O4) * (CUMcbb);
+ c_022 = -O4 * (c1o1 / omega - c1o2) * (dyuy + dzuz) * c2o3 * factorA + (c1o1 - O4) * (c_022);
+ c_202 = -O4 * (c1o1 / omega - c1o2) * (dxux + dzuz) * c2o3 * factorA + (c1o1 - O4) * (c_202);
+ c_220 = -O4 * (c1o1 / omega - c1o2) * (dyuy + dxux) * c2o3 * factorA + (c1o1 - O4) * (c_220);
+ c_112 = -O4 * (c1o1 / omega - c1o2) * Dxy * c1o3 * factorB + (c1o1 - O4) * (c_112);
+ c_121 = -O4 * (c1o1 / omega - c1o2) * Dxz * c1o3 * factorB + (c1o1 - O4) * (c_121);
+ c_211 = -O4 * (c1o1 / omega - c1o2) * Dyz * c1o3 * factorB + (c1o1 - O4) * (c_211);
+
//////////////////////////////////////////////////////////////////////////
// 5.
- CUMbcc += O5 * (-CUMbcc);
- CUMcbc += O5 * (-CUMcbc);
- CUMccb += O5 * (-CUMccb);
+ c_122 += O5 * (-c_122);
+ c_212 += O5 * (-c_212);
+ c_221 += O5 * (-c_221);
//////////////////////////////////////////////////////////////////////////
// 6.
- CUMccc += O6 * (-CUMccc);
+ c_222 += O6 * (-c_222);
////////////////////////////////////////////////////////////////////////////////////
//! - Compute central moments from post collision cumulants according to Eq. (53)-(56) in
@@ -583,68 +548,58 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//////////////////////////////////////////////////////////////////////////
// 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;
+ m_211 = c_211 + c1o3 * ((c3o1 * m_200 + c1o1) * m_011 + c6o1 * m_110 * m_101) * oneOverRho;
+ m_121 = c_121 + c1o3 * ((c3o1 * m_020 + c1o1) * m_101 + c6o1 * m_110 * m_011) * oneOverRho;
+ m_112 = c_112 + c1o3 * ((c3o1 * m_002 + c1o1) * m_110 + c6o1 * m_101 * m_011) * oneOverRho;
+
+ m_220 =
+ c_220 + (((m_200 * m_020 + c2o1 * m_110 * m_110) * c9o1 + c3o1 * (m_200 + m_020)) * oneOverRho - (drho * oneOverRho)) * c1o9;
+ m_202 =
+ c_202 + (((m_200 * m_002 + c2o1 * m_101 * m_101) * c9o1 + c3o1 * (m_200 + m_002)) * oneOverRho - (drho * oneOverRho)) * c1o9;
+ m_022 =
+ c_022 + (((m_002 * m_020 + c2o1 * m_011 * m_011) * c9o1 + c3o1 * (m_002 + m_020)) * oneOverRho - (drho * oneOverRho)) * 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;
+ m_122 = c_122 + c1o3 *
+ (c3o1 * (m_002 * m_120 + m_020 * m_102 + c4o1 * m_011 * m_111 + c2o1 * (m_101 * m_021 + m_110 * m_012)) +
+ (m_120 + m_102)) * oneOverRho;
+ m_212 = c_212 + c1o3 *
+ (c3o1 * (m_002 * m_210 + m_200 * m_012 + c4o1 * m_101 * m_111 + c2o1 * (m_011 * m_201 + m_110 * m_102)) +
+ (m_210 + m_012)) * oneOverRho;
+ m_221 = c_221 + c1o3 *
+ (c3o1 * (m_200 * m_021 + m_020 * m_201 + c4o1 * m_110 * m_111 + c2o1 * (m_101 * m_120 + m_011 * m_210)) +
+ (m_021 + m_201)) * oneOverRho;
//////////////////////////////////////////////////////////////////////////
// 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));
+ m_222 = c_222 - ((-c4o1 * m_111 * m_111 - (m_200 * m_022 + m_020 * m_202 + m_002 * m_220) -
+ c4o1 * (m_011 * m_211 + m_101 * m_121 + m_110 * m_112) -
+ c2o1 * (m_120 * m_102 + m_210 * m_012 + m_201 * m_021)) *
+ oneOverRho +
+ (c4o1 * (m_101 * m_101 * m_020 + m_011 * m_011 * m_200 + m_110 * m_110 * m_002) +
+ c2o1 * (m_200 * m_020 * m_002) + c16o1 * m_110 * m_101 * m_011) *
+ oneOverRho * oneOverRho -
+ c1o3 * (m_022 + m_202 + m_220) * oneOverRho - c1o9 * (m_200 + m_020 + m_002) * oneOverRho +
+ (c2o1 * (m_101 * m_101 + m_011 * m_011 + m_110 * m_110) +
+ (m_002 * m_020 + m_002 * m_200 + m_020 * m_200) + c1o3 * (m_002 + m_020 + m_200)) *
+ oneOverRho * oneOverRho * c2o3 +
+ c1o27 * ((drho * drho - drho) * oneOverRho * oneOverRho));
////////////////////////////////////////////////////////////////////////////////////
//! - 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;
-
+ m_100 = -m_100;
+ m_010 = -m_010;
+ m_001 = -m_001;
//Write to array here to distribute read/write
- rho[k] = drho;
- vx[k] = vvx;
- vy[k] = vvy;
- vz[k] = vvz;
+ rho[k_000] = drho;
+ vx[k_000] = vvx;
+ vy[k_000] = vvy;
+ vz[k_000] = vvz;
////////////////////////////////////////////////////////////////////////////////////
//! - Chimera transform from central moments to well conditioned distributions as defined in Appendix J in
@@ -655,39 +610,39 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//!
////////////////////////////////////////////////////////////////////////////////////
// X - Dir
- backwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1o1, c1o1);
- backwardChimera(mfaba, mfbba, mfcba, vvx, vx2);
- backwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3o1, c1o3);
- backwardChimera(mfaab, mfbab, mfcab, vvx, vx2);
- backwardChimera(mfabb, mfbbb, mfcbb, vvx, vx2);
- backwardChimera(mfacb, mfbcb, mfccb, vvx, vx2);
- backwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3o1, c1o3);
- backwardChimera(mfabc, mfbbc, mfcbc, vvx, vx2);
- backwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c9o1, c1o9);
+ backwardInverseChimeraWithK(m_000, m_100, m_200, vvx, vx2, c1o1, c1o1);
+ backwardChimera( m_010, m_110, m_210, vvx, vx2);
+ backwardInverseChimeraWithK(m_020, m_120, m_220, vvx, vx2, c3o1, c1o3);
+ backwardChimera( m_001, m_101, m_201, vvx, vx2);
+ backwardChimera( m_011, m_111, m_211, vvx, vx2);
+ backwardChimera( m_021, m_121, m_221, vvx, vx2);
+ backwardInverseChimeraWithK(m_002, m_102, m_202, vvx, vx2, c3o1, c1o3);
+ backwardChimera( m_012, m_112, m_212, vvx, vx2);
+ backwardInverseChimeraWithK(m_022, m_122, m_222, vvx, vx2, c9o1, c1o9);
////////////////////////////////////////////////////////////////////////////////////
// Y - Dir
- backwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6o1, c1o6);
- backwardChimera(mfaab, mfabb, mfacb, vvy, vy2);
- backwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18o1, c1o18);
- backwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
- backwardChimera(mfbab, mfbbb, mfbcb, vvy, vy2);
- backwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
- backwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6o1, c1o6);
- backwardChimera(mfcab, mfcbb, mfccb, vvy, vy2);
- backwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18o1, c1o18);
+ backwardInverseChimeraWithK(m_000, m_010, m_020, vvy, vy2, c6o1, c1o6);
+ backwardChimera( m_001, m_011, m_021, vvy, vy2);
+ backwardInverseChimeraWithK(m_002, m_012, m_022, vvy, vy2, c18o1, c1o18);
+ backwardInverseChimeraWithK(m_100, m_110, m_120, vvy, vy2, c3o2, c2o3);
+ backwardChimera( m_101, m_111, m_121, vvy, vy2);
+ backwardInverseChimeraWithK(m_102, m_112, m_122, vvy, vy2, c9o2, c2o9);
+ backwardInverseChimeraWithK(m_200, m_210, m_220, vvy, vy2, c6o1, c1o6);
+ backwardChimera( m_201, m_211, m_221, vvy, vy2);
+ backwardInverseChimeraWithK(m_202, m_212, m_222, vvy, vy2, c18o1, c1o18);
////////////////////////////////////////////////////////////////////////////////////
// Z - Dir
- backwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36o1, c1o36);
- backwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9o1, c1o9);
- backwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36o1, c1o36);
- backwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9o1, c1o9);
- backwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
- backwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9o1, c1o9);
- backwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36o1, c1o36);
- backwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9o1, c1o9);
- backwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36o1, c1o36);
+ backwardInverseChimeraWithK(m_000, m_001, m_002, vvz, vz2, c36o1, c1o36);
+ backwardInverseChimeraWithK(m_010, m_011, m_012, vvz, vz2, c9o1, c1o9);
+ backwardInverseChimeraWithK(m_020, m_021, m_022, vvz, vz2, c36o1, c1o36);
+ backwardInverseChimeraWithK(m_100, m_101, m_102, vvz, vz2, c9o1, c1o9);
+ backwardInverseChimeraWithK(m_110, m_111, m_112, vvz, vz2, c9o4, c4o9);
+ backwardInverseChimeraWithK(m_120, m_121, m_122, vvz, vz2, c9o1, c1o9);
+ backwardInverseChimeraWithK(m_200, m_201, m_202, vvz, vz2, c36o1, c1o36);
+ backwardInverseChimeraWithK(m_210, m_211, m_212, vvz, vz2, c9o1, c1o9);
+ backwardInverseChimeraWithK(m_220, m_221, m_222, vvz, vz2, c36o1, c1o36);
////////////////////////////////////////////////////////////////////////////////////
//! - Write distributions: style of reading and writing the distributions from/to
@@ -695,35 +650,33 @@ __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim(
//! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017),
//! DOI:10.3390/computation5020019 ]</b></a>
//!
- (dist.f[DIR_P00])[k] = mfabb;
- (dist.f[DIR_M00])[kw] = mfcbb;
- (dist.f[DIR_0P0])[k] = mfbab;
- (dist.f[DIR_0M0])[ks] = mfbcb;
- (dist.f[DIR_00P])[k] = mfbba;
- (dist.f[DIR_00M])[kb] = mfbbc;
- (dist.f[DIR_PP0])[k] = mfaab;
- (dist.f[DIR_MM0])[ksw] = mfccb;
- (dist.f[DIR_PM0])[ks] = mfacb;
- (dist.f[DIR_MP0])[kw] = mfcab;
- (dist.f[DIR_P0P])[k] = mfaba;
- (dist.f[DIR_M0M])[kbw] = mfcbc;
- (dist.f[DIR_P0M])[kb] = mfabc;
- (dist.f[DIR_M0P])[kw] = mfcba;
- (dist.f[DIR_0PP])[k] = mfbaa;
- (dist.f[DIR_0MM])[kbs] = mfbcc;
- (dist.f[DIR_0PM])[kb] = mfbac;
- (dist.f[DIR_0MP])[ks] = mfbca;
- (dist.f[DIR_000])[k] = mfbbb;
- (dist.f[DIR_PPP])[k] = mfaaa;
- (dist.f[DIR_PMP])[ks] = mfaca;
- (dist.f[DIR_PPM])[kb] = mfaac;
- (dist.f[DIR_PMM])[kbs] = mfacc;
- (dist.f[DIR_MPP])[kw] = mfcaa;
- (dist.f[DIR_MMP])[ksw] = mfcca;
- (dist.f[DIR_MPM])[kbw] = mfcac;
- (dist.f[DIR_MMM])[kbsw] = mfccc;
-
-
+ (dist.f[DIR_P00])[k_000] = f_M00;
+ (dist.f[DIR_M00])[k_M00] = f_P00;
+ (dist.f[DIR_0P0])[k_000] = f_0M0;
+ (dist.f[DIR_0M0])[k_0M0] = f_0P0;
+ (dist.f[DIR_00P])[k_000] = f_00M;
+ (dist.f[DIR_00M])[k_00M] = f_00P;
+ (dist.f[DIR_PP0])[k_000] = f_MM0;
+ (dist.f[DIR_MM0])[k_MM0] = f_PP0;
+ (dist.f[DIR_PM0])[k_0M0] = f_MP0;
+ (dist.f[DIR_MP0])[k_M00] = f_PM0;
+ (dist.f[DIR_P0P])[k_000] = f_M0M;
+ (dist.f[DIR_M0M])[k_M0M] = f_P0P;
+ (dist.f[DIR_P0M])[k_00M] = f_M0P;
+ (dist.f[DIR_M0P])[k_M00] = f_P0M;
+ (dist.f[DIR_0PP])[k_000] = f_0MM;
+ (dist.f[DIR_0MM])[k_0MM] = f_0PP;
+ (dist.f[DIR_0PM])[k_00M] = f_0MP;
+ (dist.f[DIR_0MP])[k_0M0] = f_0PM;
+ (dist.f[DIR_000])[k_000] = f_000;
+ (dist.f[DIR_PPP])[k_000] = f_MMM;
+ (dist.f[DIR_PMP])[k_0M0] = f_MPM;
+ (dist.f[DIR_PPM])[k_00M] = f_MMP;
+ (dist.f[DIR_PMM])[k_0MM] = f_MPP;
+ (dist.f[DIR_MPP])[k_M00] = f_PMM;
+ (dist.f[DIR_MMP])[k_MM0] = f_PPM;
+ (dist.f[DIR_MPM])[k_M0M] = f_PMP;
+ (dist.f[DIR_MMM])[k_MMM] = f_PPP;
}
}
@@ -731,6 +684,4 @@ template __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim < Turbu
template __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::Smagorinsky > ( real omega_in, uint* typeOfGridNode, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long size_Mat, int level, bool bodyForce, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep);
-template __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::QR > ( real omega_in, uint* typeOfGridNode, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long size_Mat, int level, bool bodyForce, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep);
-
-template __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::None > ( real omega_in, uint* typeOfGridNode, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long size_Mat, int level, bool bodyForce, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep);
\ No newline at end of file
+template __global__ void LB_Kernel_TurbulentViscosityCumulantK17CompChim < TurbulenceModel::QR > ( real omega_in, uint* typeOfGridNode, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long size_Mat, int level, bool bodyForce, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep);
\ No newline at end of file