From d688c76da91a315c6a32df1e820d28a732c187b1 Mon Sep 17 00:00:00 2001
From: Timon Habenicht <t.habenicht@tu-bs.de>
Date: Wed, 21 Nov 2018 11:04:43 +0100
Subject: [PATCH] adds Kernel CumulantF32018
---
.../CumulantF32018/CumulantF32018CompSP27.cu | 60 ++
.../CumulantF32018/CumulantF32018CompSP27.h | 22 +
.../CumulantF32018CompSP27_Device.cu | 945 ++++++++++++++++++
.../CumulantF32018CompSP27_Device.cuh | 19 +
.../CumulantF32018/package.include | 0
5 files changed, 1046 insertions(+)
create mode 100644 src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27.cu
create mode 100644 src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27.h
create mode 100644 src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27_Device.cu
create mode 100644 src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27_Device.cuh
create mode 100644 src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/package.include
diff --git a/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27.cu b/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27.cu
new file mode 100644
index 000000000..545d752f3
--- /dev/null
+++ b/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27.cu
@@ -0,0 +1,60 @@
+#include "CumulantF32018CompSP27.h"
+
+#include "CumulantF32018CompSP27_Device.cuh"
+
+#include "Parameter\Parameter.h"
+
+std::shared_ptr<Kernel> CumulantF32018CompSP27::getNewInstance(std::shared_ptr<Parameter> para, int level)
+{
+ return std::shared_ptr<Kernel>(new CumulantF32018CompSP27(para, level));
+}
+
+void CumulantF32018CompSP27::run()
+{
+ int numberOfThreads = para->getParD(level)->numberofthreads;
+ int size_Mat = para->getParD(level)->size_Mat_SP;
+
+ int Grid = (size_Mat / numberOfThreads) + 1;
+ int Grid1, Grid2;
+ if (Grid>512)
+ {
+ Grid1 = 512;
+ Grid2 = (Grid / Grid1) + 1;
+ }
+ else
+ {
+ Grid1 = 1;
+ Grid2 = Grid;
+ }
+ dim3 grid(Grid1, Grid2);
+ dim3 threads(numberOfThreads, 1, 1);
+
+ LB_Kernel_Cumulant_D3Q27F3_2018 << < grid, threads >> >( para->getParD(level)->omega,
+ para->getParD(level)->geoSP,
+ para->getParD(level)->neighborX_SP,
+ para->getParD(level)->neighborY_SP,
+ para->getParD(level)->neighborZ_SP,
+ para->getParD(level)->d0SP.f[0],
+ para->getParD(level)->g6.g[0],
+ para->getParD(level)->size_Mat_SP,
+ level,
+ para->getForcesDev(),
+ para->getParD(level)->evenOrOdd);
+ getLastCudaError("LB_Kernel_Cumulant_D3Q27F3 execution failed");
+}
+
+bool CumulantF32018CompSP27::checkParameter()
+{
+ if (para->getUseWale())
+ return false;
+ else if (para->getDiffOn())
+ return false;
+ else
+ return true;
+}
+
+CumulantF32018CompSP27::CumulantF32018CompSP27(std::shared_ptr<Parameter> para, int level)
+{
+ this->para = para;
+ this->level = level;
+}
\ No newline at end of file
diff --git a/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27.h b/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27.h
new file mode 100644
index 000000000..5712559cf
--- /dev/null
+++ b/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27.h
@@ -0,0 +1,22 @@
+#ifndef CUMULANT_F32018_COMP_SP27_H
+#define CUMULANT_F32018_COMP_SP27_H
+
+#include "Kernel\Kernel.h"
+
+#include <memory>
+
+class CumulantF32018CompSP27 : public Kernel
+{
+public:
+ static std::shared_ptr< Kernel> getNewInstance(std::shared_ptr< Parameter> para, int level);
+ void run();
+ bool checkParameter();
+
+private:
+ CumulantF32018CompSP27(std::shared_ptr< Parameter> para, int level);
+ std::shared_ptr< Parameter> para;
+ int level;
+
+};
+
+#endif
\ No newline at end of file
diff --git a/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27_Device.cu b/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27_Device.cu
new file mode 100644
index 000000000..6f5b42daa
--- /dev/null
+++ b/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27_Device.cu
@@ -0,0 +1,945 @@
+#include "CumulantF32018CompSP27_Device.cuh"
+
+#include "LBM/D3Q27.h"
+#include "math.h"
+#include "GPU/constant.h"
+
+extern "C" __global__ void LB_Kernel_Cumulant_D3Q27F3_2018(real omega,
+ unsigned int* bcMatD,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ real* DDStart,
+ real* G6,
+ int size_Mat,
+ int level,
+ real* forces,
+ bool EvenOrOdd)
+{
+ ////////////////////////////////////////////////////////////////////////////////
+ const unsigned x = threadIdx.x; // Globaler x-Index
+ const unsigned y = blockIdx.x; // Globaler y-Index
+ const unsigned z = blockIdx.y; // Globaler z-Index
+
+ const unsigned nx = blockDim.x;
+ const unsigned ny = gridDim.x;
+
+ const unsigned k = nx*(ny*z + y) + x;
+ //////////////////////////////////////////////////////////////////////////
+
+ if (k < size_Mat)
+ {
+ ////////////////////////////////////////////////////////////////////////////////
+ unsigned int BC;
+ BC = bcMatD[k];
+
+ if ((BC != GEO_SOLID) && (BC != GEO_VOID))
+ {
+ Distributions27 D;
+ if (EvenOrOdd == true)
+ {
+ D.f[dirE] = &DDStart[dirE *size_Mat];
+ D.f[dirW] = &DDStart[dirW *size_Mat];
+ D.f[dirN] = &DDStart[dirN *size_Mat];
+ D.f[dirS] = &DDStart[dirS *size_Mat];
+ D.f[dirT] = &DDStart[dirT *size_Mat];
+ D.f[dirB] = &DDStart[dirB *size_Mat];
+ D.f[dirNE] = &DDStart[dirNE *size_Mat];
+ D.f[dirSW] = &DDStart[dirSW *size_Mat];
+ D.f[dirSE] = &DDStart[dirSE *size_Mat];
+ D.f[dirNW] = &DDStart[dirNW *size_Mat];
+ D.f[dirTE] = &DDStart[dirTE *size_Mat];
+ D.f[dirBW] = &DDStart[dirBW *size_Mat];
+ D.f[dirBE] = &DDStart[dirBE *size_Mat];
+ D.f[dirTW] = &DDStart[dirTW *size_Mat];
+ D.f[dirTN] = &DDStart[dirTN *size_Mat];
+ D.f[dirBS] = &DDStart[dirBS *size_Mat];
+ D.f[dirBN] = &DDStart[dirBN *size_Mat];
+ D.f[dirTS] = &DDStart[dirTS *size_Mat];
+ D.f[dirZERO] = &DDStart[dirZERO*size_Mat];
+ D.f[dirTNE] = &DDStart[dirTNE *size_Mat];
+ D.f[dirTSW] = &DDStart[dirTSW *size_Mat];
+ D.f[dirTSE] = &DDStart[dirTSE *size_Mat];
+ D.f[dirTNW] = &DDStart[dirTNW *size_Mat];
+ D.f[dirBNE] = &DDStart[dirBNE *size_Mat];
+ D.f[dirBSW] = &DDStart[dirBSW *size_Mat];
+ D.f[dirBSE] = &DDStart[dirBSE *size_Mat];
+ D.f[dirBNW] = &DDStart[dirBNW *size_Mat];
+ }
+ else
+ {
+ D.f[dirW] = &DDStart[dirE *size_Mat];
+ D.f[dirE] = &DDStart[dirW *size_Mat];
+ D.f[dirS] = &DDStart[dirN *size_Mat];
+ D.f[dirN] = &DDStart[dirS *size_Mat];
+ D.f[dirB] = &DDStart[dirT *size_Mat];
+ D.f[dirT] = &DDStart[dirB *size_Mat];
+ D.f[dirSW] = &DDStart[dirNE *size_Mat];
+ D.f[dirNE] = &DDStart[dirSW *size_Mat];
+ D.f[dirNW] = &DDStart[dirSE *size_Mat];
+ D.f[dirSE] = &DDStart[dirNW *size_Mat];
+ D.f[dirBW] = &DDStart[dirTE *size_Mat];
+ D.f[dirTE] = &DDStart[dirBW *size_Mat];
+ D.f[dirTW] = &DDStart[dirBE *size_Mat];
+ D.f[dirBE] = &DDStart[dirTW *size_Mat];
+ D.f[dirBS] = &DDStart[dirTN *size_Mat];
+ D.f[dirTN] = &DDStart[dirBS *size_Mat];
+ D.f[dirTS] = &DDStart[dirBN *size_Mat];
+ D.f[dirBN] = &DDStart[dirTS *size_Mat];
+ D.f[dirZERO] = &DDStart[dirZERO*size_Mat];
+ D.f[dirBSW] = &DDStart[dirTNE *size_Mat];
+ D.f[dirBNE] = &DDStart[dirTSW *size_Mat];
+ D.f[dirBNW] = &DDStart[dirTSE *size_Mat];
+ D.f[dirBSE] = &DDStart[dirTNW *size_Mat];
+ D.f[dirTSW] = &DDStart[dirBNE *size_Mat];
+ D.f[dirTNE] = &DDStart[dirBSW *size_Mat];
+ D.f[dirTNW] = &DDStart[dirBSE *size_Mat];
+ D.f[dirTSE] = &DDStart[dirBNW *size_Mat];
+ }
+
+ Distributions6 G;
+ if (EvenOrOdd == true)
+ {
+ G.g[dirE] = &G6[dirE *size_Mat];
+ G.g[dirW] = &G6[dirW *size_Mat];
+ G.g[dirN] = &G6[dirN *size_Mat];
+ G.g[dirS] = &G6[dirS *size_Mat];
+ G.g[dirT] = &G6[dirT *size_Mat];
+ G.g[dirB] = &G6[dirB *size_Mat];
+ }
+ else
+ {
+ G.g[dirW] = &G6[dirE *size_Mat];
+ G.g[dirE] = &G6[dirW *size_Mat];
+ G.g[dirS] = &G6[dirN *size_Mat];
+ G.g[dirN] = &G6[dirS *size_Mat];
+ G.g[dirB] = &G6[dirT *size_Mat];
+ G.g[dirT] = &G6[dirB *size_Mat];
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////
+ //index
+ //unsigned int kzero= k;
+ //unsigned int ke = k;
+ unsigned int kw = neighborX[k];
+ //unsigned int kn = k;
+ unsigned int ks = neighborY[k];
+ //unsigned int kt = k;
+ unsigned int kb = neighborZ[k];
+ unsigned int ksw = neighborY[kw];
+ //unsigned int kne = k;
+ //unsigned int kse = ks;
+ //unsigned int knw = kw;
+ unsigned int kbw = neighborZ[kw];
+ //unsigned int kte = k;
+ //unsigned int kbe = kb;
+ //unsigned int ktw = kw;
+ unsigned int kbs = neighborZ[ks];
+ //unsigned int ktn = k;
+ //unsigned int kbn = kb;
+ //unsigned int kts = ks;
+ //unsigned int ktse = ks;
+ //unsigned int kbnw = kbw;
+ //unsigned int ktnw = kw;
+ //unsigned int kbse = kbs;
+ //unsigned int ktsw = ksw;
+ //unsigned int kbne = kb;
+ //unsigned int ktne = k;
+ unsigned int kbsw = neighborZ[ksw];
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ real mgcbb = (G.g[dirE])[k];
+ real mgabb = (G.g[dirW])[kw];
+ real mgbcb = (G.g[dirN])[k];
+ real mgbab = (G.g[dirS])[ks];
+ real mgbbc = (G.g[dirT])[k];
+ real mgbba = (G.g[dirB])[kb];
+ real dxuxdxux = c1o2 * (-mgcbb + mgabb);
+ real dyuydyuy = c1o2 * (-mgbcb + mgbab);
+ real dzuzdzuz = c1o2 * (-mgbbc + mgbba);
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ real mfcbb = (D.f[dirE])[k];
+ real mfabb = (D.f[dirW])[kw];
+ real mfbcb = (D.f[dirN])[k];
+ real mfbab = (D.f[dirS])[ks];
+ real mfbbc = (D.f[dirT])[k];
+ real mfbba = (D.f[dirB])[kb];
+ real mfccb = (D.f[dirNE])[k];
+ real mfaab = (D.f[dirSW])[ksw];
+ real mfcab = (D.f[dirSE])[ks];
+ real mfacb = (D.f[dirNW])[kw];
+ real mfcbc = (D.f[dirTE])[k];
+ real mfaba = (D.f[dirBW])[kbw];
+ real mfcba = (D.f[dirBE])[kb];
+ real mfabc = (D.f[dirTW])[kw];
+ real mfbcc = (D.f[dirTN])[k];
+ real mfbaa = (D.f[dirBS])[kbs];
+ real mfbca = (D.f[dirBN])[kb];
+ real mfbac = (D.f[dirTS])[ks];
+ real mfbbb = (D.f[dirZERO])[k];
+ real mfccc = (D.f[dirTNE])[k];
+ real mfaac = (D.f[dirTSW])[ksw];
+ real mfcac = (D.f[dirTSE])[ks];
+ real mfacc = (D.f[dirTNW])[kw];
+ real mfcca = (D.f[dirBNE])[kb];
+ real mfaaa = (D.f[dirBSW])[kbsw];
+ real mfcaa = (D.f[dirBSE])[kbs];
+ real mfaca = (D.f[dirBNW])[kbw];
+ ////////////////////////////////////////////////////////////////////////////////////
+ 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 = one + drho;
+ ////////////////////////////////////////////////////////////////////////////////////
+ real vvx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+ (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
+ (mfcbb - mfabb)) / rho;
+ real vvy = ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+ (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
+ (mfbcb - mfbab)) / rho;
+ real vvz = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+ (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
+ (mfbbc - mfbba)) / rho;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //the force be with you
+ real fx = forces[0] / (pow(two, level)); //zero;//0.0032653/(pow(two,level)); //0.000000005;//(two/1600000.0) / 120.0; //
+ real fy = forces[1] / (pow(two, level)); //zero;
+ real fz = forces[2] / (pow(two, level)); //zero;
+ vvx += fx;
+ vvy += fy;
+ vvz += fz;
+ ////////////////////////////////////////////////////////////////////////////////////
+ real oMdrho = one; // comp special
+ real m0, m1, m2;
+ real vx2;
+ real vy2;
+ real vz2;
+ vx2 = vvx*vvx;
+ vy2 = vvy*vvy;
+ vz2 = vvz*vvz;
+ ////////////////////////////////////////////////////////////////////////////////////
+ real wadjust;
+ real qudricLimitP = 0.01f;// * 0.0001f;
+ real qudricLimitM = 0.01f;// * 0.0001f;
+ real qudricLimitD = 0.01f;// * 0.001f;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //Hin
+ ////////////////////////////////////////////////////////////////////////////////////
+ // mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ m2 = mfaaa + mfaac;
+ m1 = mfaac - mfaaa;
+ m0 = m2 + mfaab;
+ mfaaa = m0;
+ m0 += c1o36 * oMdrho;
+ mfaab = m1 - m0 * vvz;
+ mfaac = m2 - two* m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaba + mfabc;
+ m1 = mfabc - mfaba;
+ m0 = m2 + mfabb;
+ mfaba = m0;
+ m0 += c1o9 * oMdrho;
+ mfabb = m1 - m0 * vvz;
+ mfabc = m2 - two* m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaca + mfacc;
+ m1 = mfacc - mfaca;
+ m0 = m2 + mfacb;
+ mfaca = m0;
+ m0 += c1o36 * oMdrho;
+ mfacb = m1 - m0 * vvz;
+ mfacc = m2 - two* m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbaa + mfbac;
+ m1 = mfbac - mfbaa;
+ m0 = m2 + mfbab;
+ mfbaa = m0;
+ m0 += c1o9 * oMdrho;
+ mfbab = m1 - m0 * vvz;
+ mfbac = m2 - two* m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbba + mfbbc;
+ m1 = mfbbc - mfbba;
+ m0 = m2 + mfbbb;
+ mfbba = m0;
+ m0 += c4o9 * oMdrho;
+ mfbbb = m1 - m0 * vvz;
+ mfbbc = m2 - two* m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbca + mfbcc;
+ m1 = mfbcc - mfbca;
+ m0 = m2 + mfbcb;
+ mfbca = m0;
+ m0 += c1o9 * oMdrho;
+ mfbcb = m1 - m0 * vvz;
+ mfbcc = m2 - two* m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcaa + mfcac;
+ m1 = mfcac - mfcaa;
+ m0 = m2 + mfcab;
+ mfcaa = m0;
+ m0 += c1o36 * oMdrho;
+ mfcab = m1 - m0 * vvz;
+ mfcac = m2 - two* m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcba + mfcbc;
+ m1 = mfcbc - mfcba;
+ m0 = m2 + mfcbb;
+ mfcba = m0;
+ m0 += c1o9 * oMdrho;
+ mfcbb = m1 - m0 * vvz;
+ mfcbc = m2 - two* m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcca + mfccc;
+ m1 = mfccc - mfcca;
+ m0 = m2 + mfccb;
+ mfcca = m0;
+ m0 += c1o36 * oMdrho;
+ mfccb = m1 - m0 * vvz;
+ mfccc = m2 - two* m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ // mit 1/6, 0, 1/18, 2/3, 0, 2/9, 1/6, 0, 1/18 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ m2 = mfaaa + mfaca;
+ m1 = mfaca - mfaaa;
+ m0 = m2 + mfaba;
+ mfaaa = m0;
+ m0 += c1o6 * oMdrho;
+ mfaba = m1 - m0 * vvy;
+ mfaca = m2 - two* m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaab + mfacb;
+ m1 = mfacb - mfaab;
+ m0 = m2 + mfabb;
+ mfaab = m0;
+ mfabb = m1 - m0 * vvy;
+ mfacb = m2 - two* m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaac + mfacc;
+ m1 = mfacc - mfaac;
+ m0 = m2 + mfabc;
+ mfaac = m0;
+ m0 += c1o18 * oMdrho;
+ mfabc = m1 - m0 * vvy;
+ mfacc = m2 - two* m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbaa + mfbca;
+ m1 = mfbca - mfbaa;
+ m0 = m2 + mfbba;
+ mfbaa = m0;
+ m0 += c2o3 * oMdrho;
+ mfbba = m1 - m0 * vvy;
+ mfbca = m2 - two* m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbab + mfbcb;
+ m1 = mfbcb - mfbab;
+ m0 = m2 + mfbbb;
+ mfbab = m0;
+ mfbbb = m1 - m0 * vvy;
+ mfbcb = m2 - two* m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbac + mfbcc;
+ m1 = mfbcc - mfbac;
+ m0 = m2 + mfbbc;
+ mfbac = m0;
+ m0 += c2o9 * oMdrho;
+ mfbbc = m1 - m0 * vvy;
+ mfbcc = m2 - two* m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcaa + mfcca;
+ m1 = mfcca - mfcaa;
+ m0 = m2 + mfcba;
+ mfcaa = m0;
+ m0 += c1o6 * oMdrho;
+ mfcba = m1 - m0 * vvy;
+ mfcca = m2 - two* m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcab + mfccb;
+ m1 = mfccb - mfcab;
+ m0 = m2 + mfcbb;
+ mfcab = m0;
+ mfcbb = m1 - m0 * vvy;
+ mfccb = m2 - two* m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcac + mfccc;
+ m1 = mfccc - mfcac;
+ m0 = m2 + mfcbc;
+ mfcac = m0;
+ m0 += c1o18 * oMdrho;
+ mfcbc = m1 - m0 * vvy;
+ mfccc = m2 - two* m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ // mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ m2 = mfaaa + mfcaa;
+ m1 = mfcaa - mfaaa;
+ m0 = m2 + mfbaa;
+ mfaaa = m0;
+ m0 += one* oMdrho;
+ mfbaa = m1 - m0 * vvx;
+ mfcaa = m2 - two* m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaba + mfcba;
+ m1 = mfcba - mfaba;
+ m0 = m2 + mfbba;
+ mfaba = m0;
+ mfbba = m1 - m0 * vvx;
+ mfcba = m2 - two* m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaca + mfcca;
+ m1 = mfcca - mfaca;
+ m0 = m2 + mfbca;
+ mfaca = m0;
+ m0 += c1o3 * oMdrho;
+ mfbca = m1 - m0 * vvx;
+ mfcca = m2 - two* m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaab + mfcab;
+ m1 = mfcab - mfaab;
+ m0 = m2 + mfbab;
+ mfaab = m0;
+ mfbab = m1 - m0 * vvx;
+ mfcab = m2 - two* m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfabb + mfcbb;
+ m1 = mfcbb - mfabb;
+ m0 = m2 + mfbbb;
+ mfabb = m0;
+ mfbbb = m1 - m0 * vvx;
+ mfcbb = m2 - two* m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfacb + mfccb;
+ m1 = mfccb - mfacb;
+ m0 = m2 + mfbcb;
+ mfacb = m0;
+ mfbcb = m1 - m0 * vvx;
+ mfccb = m2 - two* m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaac + mfcac;
+ m1 = mfcac - mfaac;
+ m0 = m2 + mfbac;
+ mfaac = m0;
+ m0 += c1o3 * oMdrho;
+ mfbac = m1 - m0 * vvx;
+ mfcac = m2 - two* m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfabc + mfcbc;
+ m1 = mfcbc - mfabc;
+ m0 = m2 + mfbbc;
+ mfabc = m0;
+ mfbbc = m1 - m0 * vvx;
+ mfcbc = m2 - two* m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfacc + mfccc;
+ m1 = mfccc - mfacc;
+ m0 = m2 + mfbcc;
+ mfacc = m0;
+ m0 += c1o9 * oMdrho;
+ mfbcc = m1 - m0 * vvx;
+ mfccc = m2 - two* m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Cumulants
+ ////////////////////////////////////////////////////////////////////////////////////
+ real OxxPyyPzz = one; //set the bulk viscosity one is high / two is very low and zero is (too) high ... (also called omega 2)
+
+ ////////////////////////////////////////////////////////////
+ //3.
+ //////////////////////////////
+ real OxyyPxzz = eight*(-two + omega)*(one + two*omega) / (-eight - fourteen*omega + seven*omega*omega);//one;
+ real OxyyMxzz = eight*(-two + omega)*(-seven + four*omega) / (fiftysix - fifty*omega + nine*omega*omega);//one;
+ real Oxyz = twentyfour*(-two + omega)*(-two - seven*omega + three*omega*omega) / (fourtyeight + c152*omega - c130*omega*omega + twentynine*omega*omega*omega);//one;
+ ////////////////////////////////////////////////////////////
+ //4.
+ //////////////////////////////
+ real O4 = one;
+ //////////////////////////////
+ //real O4 = omega;//TRT
+ ////////////////////////////////////////////////////////////
+ //5.
+ //////////////////////////////
+ real O5 = one;
+ ////////////////////////////////////////////////////////////
+ //6.
+ //////////////////////////////
+ real O6 = one;
+ ////////////////////////////////////////////////////////////
+
+
+ //central moments to cumulants
+ //4.
+ real CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + two * mfbba * mfbab) / rho;
+ real CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + two * mfbba * mfabb) / rho;
+ real CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + two * mfbab * mfabb) / rho;
+
+ real CUMcca = mfcca - (((mfcaa * mfaca + two * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho - c1o9*(drho / rho));
+ real CUMcac = mfcac - (((mfcaa * mfaac + two * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho - c1o9*(drho / rho));
+ real CUMacc = mfacc - (((mfaac * mfaca + two * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho - c1o9*(drho / rho));
+
+ //5.
+ real CUMbcc = mfbcc - ((mfaac * mfbca + mfaca * mfbac + four * mfabb * mfbbb + two * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac)) / rho;
+ real CUMcbc = mfcbc - ((mfaac * mfcba + mfcaa * mfabc + four * mfbab * mfbbb + two * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc)) / rho;
+ real CUMccb = mfccb - ((mfcaa * mfacb + mfaca * mfcab + four * mfbba * mfbbb + two * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab)) / rho;
+
+ //6.
+
+ real CUMccc = mfccc + ((-four * mfbbb * mfbbb
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - four * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+ - two * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) / rho
+ + (four * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + two * (mfcaa * mfaca * mfaac)
+ + sixteen * mfbba * mfbab * mfabb) / (rho * rho)
+ - c1o3 * (mfacc + mfcac + mfcca) / rho
+ - c1o9 * (mfcaa + mfaca + mfaac) / rho
+ + (two * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
+ + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) / (rho * rho) * c2o3
+ + c1o27*((drho * drho - drho) / (rho*rho)));
+
+ //2.
+ // linear combinations
+ real mxxPyyPzz = mfcaa + mfaca + mfaac;
+ real mxxMyy = mfcaa - mfaca;
+ real mxxMzz = mfcaa - mfaac;
+
+ ////////////////////////////////////////////////////////////////////////////
+ real Dxy = -three*omega*mfbba;
+ real Dxz = -three*omega*mfbab;
+ real Dyz = -three*omega*mfabb;
+
+ //3.
+ // linear combinations
+
+ 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 (hat noch nicht so gut funktioniert...Optimierungsbedarf??)
+
+ real dxux = c1o2 * (-omega) *(mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
+ real dyuy = dxux + omega * c3o2 * mxxMyy;
+ real dzuz = dxux + omega * c3o2 * mxxMzz;
+
+ mgabb = vvx*dxux;
+ mgbab = vvy*dyuy;
+ mgbba = vvz*dzuz;
+
+ mgcbb = vvx*dxux;
+ mgbcb = vvy*dyuy;
+ mgbbc = vvz*dzuz;
+
+ //relax
+ mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz) - three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz)
+ + (six - three * (omega + OxxPyyPzz) + omega * OxxPyyPzz) / (three * omega) *
+ (dxuxdxux + dyuydyuy + dzuzdzuz);
+ mxxMyy += omega * (-mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy)
+ + omega * (two*(one / omega - c1o2) * (one / omega - c1o2) - c1o6) * (dxuxdxux - dyuydyuy);
+ mxxMzz += omega * (-mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz)
+ + omega * (two*(one / omega - c1o2) * (one / omega - c1o2) - c1o6) *(dxuxdxux - dzuzdzuz);
+
+ ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////no correction
+ //mxxPyyPzz += OxxPyyPzz*(mfaaa-mxxPyyPzz);//-magicBulk*OxxPyyPzz;
+ //mxxMyy += -(-omega) * (-mxxMyy);
+ //mxxMzz += -(-omega) * (-mxxMzz);
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ mfabb += omega * (-mfabb);
+ mfbab += omega * (-mfbab);
+ mfbba += omega * (-mfbba);
+ //////////////////////////////////////////////////////////////////////////
+
+ // linear combinations back
+ mfcaa = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz);
+ mfaca = c1o3 * (-two* mxxMyy + mxxMzz + mxxPyyPzz);
+ mfaac = c1o3 * (mxxMyy - two* mxxMzz + mxxPyyPzz);
+
+
+ //relax
+ //////////////////////////////////////////////////////////////////////////
+ //das ist der limiter
+ wadjust = Oxyz + (one - Oxyz)*abs(mfbbb) / (abs(mfbbb) + qudricLimitD);
+ mfbbb += wadjust * (-mfbbb);
+ wadjust = OxyyPxzz + (one - OxyyPxzz)*abs(mxxyPyzz) / (abs(mxxyPyzz) + qudricLimitP);
+ mxxyPyzz += wadjust * (-mxxyPyzz);
+ wadjust = OxyyMxzz + (one - OxyyMxzz)*abs(mxxyMyzz) / (abs(mxxyMyzz) + qudricLimitM);
+ mxxyMyzz += wadjust * (-mxxyMyzz);
+ wadjust = OxyyPxzz + (one - OxyyPxzz)*abs(mxxzPyyz) / (abs(mxxzPyyz) + qudricLimitP);
+ mxxzPyyz += wadjust * (-mxxzPyyz);
+ wadjust = OxyyMxzz + (one - OxyyMxzz)*abs(mxxzMyyz) / (abs(mxxzMyyz) + qudricLimitM);
+ mxxzMyyz += wadjust * (-mxxzMyyz);
+ wadjust = OxyyPxzz + (one - OxyyPxzz)*abs(mxyyPxzz) / (abs(mxyyPxzz) + qudricLimitP);
+ mxyyPxzz += wadjust * (-mxyyPxzz);
+ wadjust = OxyyMxzz + (one - OxyyMxzz)*abs(mxyyMxzz) / (abs(mxyyMxzz) + qudricLimitM);
+ mxyyMxzz += wadjust * (-mxyyMxzz);
+ //////////////////////////////////////////////////////////////////////////
+ //ohne limiter
+ //mfbbb += OxyyMxzz * (-mfbbb);
+ //mxxyPyzz += OxyyPxzz * (-mxxyPyzz);
+ //mxxyMyzz += OxyyMxzz * (-mxxyMyzz);
+ //mxxzPyyz += OxyyPxzz * (-mxxzPyyz);
+ //mxxzMyyz += OxyyMxzz * (-mxxzMyyz);
+ //mxyyPxzz += OxyyPxzz * (-mxyyPxzz);
+ //mxyyMxzz += OxyyMxzz * (-mxyyMxzz);
+ //////////////////////////////////////////////////////////////////////////
+
+ // linear combinations back
+ 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.
+ //////////////////////////////////////////////////////////////////////////
+ //mit limiter
+ // wadjust = O4+(one-O4)*abs(CUMacc)/(abs(CUMacc)+qudricLimit);
+ //CUMacc += wadjust * (-CUMacc);
+ // wadjust = O4+(one-O4)*abs(CUMcac)/(abs(CUMcac)+qudricLimit);
+ //CUMcac += wadjust * (-CUMcac);
+ // wadjust = O4+(one-O4)*abs(CUMcca)/(abs(CUMcca)+qudricLimit);
+ //CUMcca += wadjust * (-CUMcca);
+
+ // wadjust = O4+(one-O4)*abs(CUMbbc)/(abs(CUMbbc)+qudricLimit);
+ //CUMbbc += wadjust * (-CUMbbc);
+ // wadjust = O4+(one-O4)*abs(CUMbcb)/(abs(CUMbcb)+qudricLimit);
+ //CUMbcb += wadjust * (-CUMbcb);
+ // wadjust = O4+(one-O4)*abs(CUMcbb)/(abs(CUMcbb)+qudricLimit);
+ //CUMcbb += wadjust * (-CUMcbb);
+ //////////////////////////////////////////////////////////////////////////
+ //ohne limiter
+ //CUMacc += O4 * (-CUMacc);
+ //CUMcac += O4 * (-CUMcac);
+ //CUMcca += O4 * (-CUMcca);
+ //CUMbbc += O4 * (-CUMbbc);
+ //CUMbcb += O4 * (-CUMbcb);
+ //CUMcbb += O4 * (-CUMcbb);
+ CUMacc = -O4*(one / omega - c1o2)*(dyuy + dzuz)*c2o3 *(four + two*omega - three*omega*omega) / (two - seven*omega + five*omega*omega) + (one - O4) * (CUMacc);
+ CUMcac = -O4*(one / omega - c1o2)*(dxux + dzuz)*c2o3 *(four + two*omega - three*omega*omega) / (two - seven*omega + five*omega*omega) + (one - O4) * (CUMcac);
+ CUMcca = -O4*(one / omega - c1o2)*(dyuy + dxux)*c2o3 *(four + two*omega - three*omega*omega) / (two - seven*omega + five*omega*omega) + (one - O4) * (CUMcca);
+ CUMbbc = -O4*(one / omega - c1o2)*Dxy*c1o3 *(four + twentyeight*omega - fourteen*omega*omega) / (six - twentyone*omega + fiveteen*omega*omega) + (one - O4) * (CUMbbc);
+ CUMbcb = -O4*(one / omega - c1o2)*Dxz*c1o3 *(four + twentyeight*omega - fourteen*omega*omega) / (six - twentyone*omega + fiveteen*omega*omega) + (one - O4) * (CUMbcb);
+ CUMcbb = -O4*(one / omega - c1o2)*Dyz*c1o3 *(four + twentyeight*omega - fourteen*omega*omega) / (six - twentyone*omega + fiveteen*omega*omega) + (one - O4) * (CUMcbb);
+ //////////////////////////////////////////////////////////////////////////
+
+
+ //5.
+ CUMbcc += O5 * (-CUMbcc);
+ CUMcbc += O5 * (-CUMcbc);
+ CUMccb += O5 * (-CUMccb);
+
+ //6.
+ CUMccc += O6 * (-CUMccc);
+
+
+
+ //back cumulants to central moments
+ //4.
+ mfcbb = CUMcbb + ((mfcaa + c1o3) * mfabb + two * mfbba * mfbab) / rho;
+ mfbcb = CUMbcb + ((mfaca + c1o3) * mfbab + two * mfbba * mfabb) / rho;
+ mfbbc = CUMbbc + ((mfaac + c1o3) * mfbba + two * mfbab * mfabb) / rho;
+
+ mfcca = CUMcca + (((mfcaa * mfaca + two * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho - c1o9*(drho / rho));
+ mfcac = CUMcac + (((mfcaa * mfaac + two * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho - c1o9*(drho / rho));
+ mfacc = CUMacc + (((mfaac * mfaca + two * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho - c1o9*(drho / rho));
+
+ //5.
+ mfbcc = CUMbcc + ((mfaac * mfbca + mfaca * mfbac + four * mfabb * mfbbb + two * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac)) / rho;
+ mfcbc = CUMcbc + ((mfaac * mfcba + mfcaa * mfabc + four * mfbab * mfbbb + two * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc)) / rho;
+ mfccb = CUMccb + ((mfcaa * mfacb + mfaca * mfcab + four * mfbba * mfbbb + two * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab)) / rho;
+
+ //6.
+
+ mfccc = CUMccc - ((-four * mfbbb * mfbbb
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - four * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+ - two * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) / rho
+ + (four * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + two * (mfcaa * mfaca * mfaac)
+ + sixteen * mfbba * mfbab * mfabb) / (rho * rho)
+ - c1o3 * (mfacc + mfcac + mfcca) / rho
+ - c1o9 * (mfcaa + mfaca + mfaac) / rho
+ + (two * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
+ + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) / (rho * rho) * c2o3
+ + c1o27*((drho * drho - drho) / (rho*rho)));
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //the force be with you
+ mfbaa = -mfbaa;
+ mfaba = -mfaba;
+ mfaab = -mfaab;
+ ////////////////////////////////////////////////////////////////////////////////////
+
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //back
+ ////////////////////////////////////////////////////////////////////////////////////
+ //mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ m0 = mfaac * c1o2 + mfaab * (vvz - c1o2) + (mfaaa + one* oMdrho) * (vz2 - vvz) * c1o2;
+ m1 = -mfaac - two* mfaab * vvz + mfaaa * (one - vz2) - one* oMdrho * vz2;
+ m2 = mfaac * c1o2 + mfaab * (vvz + c1o2) + (mfaaa + one* oMdrho) * (vz2 + vvz) * c1o2;
+ mfaaa = m0;
+ mfaab = m1;
+ mfaac = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfabc * c1o2 + mfabb * (vvz - c1o2) + mfaba * (vz2 - vvz) * c1o2;
+ m1 = -mfabc - two* mfabb * vvz + mfaba * (one - vz2);
+ m2 = mfabc * c1o2 + mfabb * (vvz + c1o2) + mfaba * (vz2 + vvz) * c1o2;
+ mfaba = m0;
+ mfabb = m1;
+ mfabc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfacc * c1o2 + mfacb * (vvz - c1o2) + (mfaca + c1o3 * oMdrho) * (vz2 - vvz) * c1o2;
+ m1 = -mfacc - two* mfacb * vvz + mfaca * (one - vz2) - c1o3 * oMdrho * vz2;
+ m2 = mfacc * c1o2 + mfacb * (vvz + c1o2) + (mfaca + c1o3 * oMdrho) * (vz2 + vvz) * c1o2;
+ mfaca = m0;
+ mfacb = m1;
+ mfacc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfbac * c1o2 + mfbab * (vvz - c1o2) + mfbaa * (vz2 - vvz) * c1o2;
+ m1 = -mfbac - two* mfbab * vvz + mfbaa * (one - vz2);
+ m2 = mfbac * c1o2 + mfbab * (vvz + c1o2) + mfbaa * (vz2 + vvz) * c1o2;
+ mfbaa = m0;
+ mfbab = m1;
+ mfbac = m2;
+ /////////b//////////////////////////////////////////////////////////////////////////
+ m0 = mfbbc * c1o2 + mfbbb * (vvz - c1o2) + mfbba * (vz2 - vvz) * c1o2;
+ m1 = -mfbbc - two* mfbbb * vvz + mfbba * (one - vz2);
+ m2 = mfbbc * c1o2 + mfbbb * (vvz + c1o2) + mfbba * (vz2 + vvz) * c1o2;
+ mfbba = m0;
+ mfbbb = m1;
+ mfbbc = m2;
+ /////////b//////////////////////////////////////////////////////////////////////////
+ m0 = mfbcc * c1o2 + mfbcb * (vvz - c1o2) + mfbca * (vz2 - vvz) * c1o2;
+ m1 = -mfbcc - two* mfbcb * vvz + mfbca * (one - vz2);
+ m2 = mfbcc * c1o2 + mfbcb * (vvz + c1o2) + mfbca * (vz2 + vvz) * c1o2;
+ mfbca = m0;
+ mfbcb = m1;
+ mfbcc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcac * c1o2 + mfcab * (vvz - c1o2) + (mfcaa + c1o3 * oMdrho) * (vz2 - vvz) * c1o2;
+ m1 = -mfcac - two* mfcab * vvz + mfcaa * (one - vz2) - c1o3 * oMdrho * vz2;
+ m2 = mfcac * c1o2 + mfcab * (vvz + c1o2) + (mfcaa + c1o3 * oMdrho) * (vz2 + vvz) * c1o2;
+ mfcaa = m0;
+ mfcab = m1;
+ mfcac = m2;
+ /////////c//////////////////////////////////////////////////////////////////////////
+ m0 = mfcbc * c1o2 + mfcbb * (vvz - c1o2) + mfcba * (vz2 - vvz) * c1o2;
+ m1 = -mfcbc - two* mfcbb * vvz + mfcba * (one - vz2);
+ m2 = mfcbc * c1o2 + mfcbb * (vvz + c1o2) + mfcba * (vz2 + vvz) * c1o2;
+ mfcba = m0;
+ mfcbb = m1;
+ mfcbc = m2;
+ /////////c//////////////////////////////////////////////////////////////////////////
+ m0 = mfccc * c1o2 + mfccb * (vvz - c1o2) + (mfcca + c1o9 * oMdrho) * (vz2 - vvz) * c1o2;
+ m1 = -mfccc - two* mfccb * vvz + mfcca * (one - vz2) - c1o9 * oMdrho * vz2;
+ m2 = mfccc * c1o2 + mfccb * (vvz + c1o2) + (mfcca + c1o9 * oMdrho) * (vz2 + vvz) * c1o2;
+ mfcca = m0;
+ mfccb = m1;
+ mfccc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ //mit 1/6, 2/3, 1/6, 0, 0, 0, 1/18, 2/9, 1/18 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ m0 = mfaca * c1o2 + mfaba * (vvy - c1o2) + (mfaaa + c1o6 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfaca - two* mfaba * vvy + mfaaa * (one - vy2) - c1o6 * oMdrho * vy2;
+ m2 = mfaca * c1o2 + mfaba * (vvy + c1o2) + (mfaaa + c1o6 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfaaa = m0;
+ mfaba = m1;
+ mfaca = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfacb * c1o2 + mfabb * (vvy - c1o2) + (mfaab + c2o3 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfacb - two* mfabb * vvy + mfaab * (one - vy2) - c2o3 * oMdrho * vy2;
+ m2 = mfacb * c1o2 + mfabb * (vvy + c1o2) + (mfaab + c2o3 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfaab = m0;
+ mfabb = m1;
+ mfacb = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfacc * c1o2 + mfabc * (vvy - c1o2) + (mfaac + c1o6 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfacc - two* mfabc * vvy + mfaac * (one - vy2) - c1o6 * oMdrho * vy2;
+ m2 = mfacc * c1o2 + mfabc * (vvy + c1o2) + (mfaac + c1o6 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfaac = m0;
+ mfabc = m1;
+ mfacc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfbca * c1o2 + mfbba * (vvy - c1o2) + mfbaa * (vy2 - vvy) * c1o2;
+ m1 = -mfbca - two* mfbba * vvy + mfbaa * (one - vy2);
+ m2 = mfbca * c1o2 + mfbba * (vvy + c1o2) + mfbaa * (vy2 + vvy) * c1o2;
+ mfbaa = m0;
+ mfbba = m1;
+ mfbca = m2;
+ /////////b//////////////////////////////////////////////////////////////////////////
+ m0 = mfbcb * c1o2 + mfbbb * (vvy - c1o2) + mfbab * (vy2 - vvy) * c1o2;
+ m1 = -mfbcb - two* mfbbb * vvy + mfbab * (one - vy2);
+ m2 = mfbcb * c1o2 + mfbbb * (vvy + c1o2) + mfbab * (vy2 + vvy) * c1o2;
+ mfbab = m0;
+ mfbbb = m1;
+ mfbcb = m2;
+ /////////b//////////////////////////////////////////////////////////////////////////
+ m0 = mfbcc * c1o2 + mfbbc * (vvy - c1o2) + mfbac * (vy2 - vvy) * c1o2;
+ m1 = -mfbcc - two* mfbbc * vvy + mfbac * (one - vy2);
+ m2 = mfbcc * c1o2 + mfbbc * (vvy + c1o2) + mfbac * (vy2 + vvy) * c1o2;
+ mfbac = m0;
+ mfbbc = m1;
+ mfbcc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcca * c1o2 + mfcba * (vvy - c1o2) + (mfcaa + c1o18 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfcca - two* mfcba * vvy + mfcaa * (one - vy2) - c1o18 * oMdrho * vy2;
+ m2 = mfcca * c1o2 + mfcba * (vvy + c1o2) + (mfcaa + c1o18 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfcaa = m0;
+ mfcba = m1;
+ mfcca = m2;
+ /////////c//////////////////////////////////////////////////////////////////////////
+ m0 = mfccb * c1o2 + mfcbb * (vvy - c1o2) + (mfcab + c2o9 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfccb - two* mfcbb * vvy + mfcab * (one - vy2) - c2o9 * oMdrho * vy2;
+ m2 = mfccb * c1o2 + mfcbb * (vvy + c1o2) + (mfcab + c2o9 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfcab = m0;
+ mfcbb = m1;
+ mfccb = m2;
+ /////////c//////////////////////////////////////////////////////////////////////////
+ m0 = mfccc * c1o2 + mfcbc * (vvy - c1o2) + (mfcac + c1o18 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfccc - two* mfcbc * vvy + mfcac * (one - vy2) - c1o18 * oMdrho * vy2;
+ m2 = mfccc * c1o2 + mfcbc * (vvy + c1o2) + (mfcac + c1o18 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfcac = m0;
+ mfcbc = m1;
+ mfccc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ //mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ m0 = mfcaa * c1o2 + mfbaa * (vvx - c1o2) + (mfaaa + c1o36 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcaa - two* mfbaa * vvx + mfaaa * (one - vx2) - c1o36 * oMdrho * vx2;
+ m2 = mfcaa * c1o2 + mfbaa * (vvx + c1o2) + (mfaaa + c1o36 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaaa = m0;
+ mfbaa = m1;
+ mfcaa = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcba * c1o2 + mfbba * (vvx - c1o2) + (mfaba + c1o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcba - two* mfbba * vvx + mfaba * (one - vx2) - c1o9 * oMdrho * vx2;
+ m2 = mfcba * c1o2 + mfbba * (vvx + c1o2) + (mfaba + c1o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaba = m0;
+ mfbba = m1;
+ mfcba = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcca * c1o2 + mfbca * (vvx - c1o2) + (mfaca + c1o36 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcca - two* mfbca * vvx + mfaca * (one - vx2) - c1o36 * oMdrho * vx2;
+ m2 = mfcca * c1o2 + mfbca * (vvx + c1o2) + (mfaca + c1o36 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaca = m0;
+ mfbca = m1;
+ mfcca = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcab * c1o2 + mfbab * (vvx - c1o2) + (mfaab + c1o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcab - two* mfbab * vvx + mfaab * (one - vx2) - c1o9 * oMdrho * vx2;
+ m2 = mfcab * c1o2 + mfbab * (vvx + c1o2) + (mfaab + c1o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaab = m0;
+ mfbab = m1;
+ mfcab = m2;
+ ///////////b////////////////////////////////////////////////////////////////////////
+ m0 = mfcbb * c1o2 + mfbbb * (vvx - c1o2) + (mfabb + c4o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcbb - two* mfbbb * vvx + mfabb * (one - vx2) - c4o9 * oMdrho * vx2;
+ m2 = mfcbb * c1o2 + mfbbb * (vvx + c1o2) + (mfabb + c4o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfabb = m0;
+ mfbbb = m1;
+ mfcbb = m2;
+ ///////////b////////////////////////////////////////////////////////////////////////
+ m0 = mfccb * c1o2 + mfbcb * (vvx - c1o2) + (mfacb + c1o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfccb - two* mfbcb * vvx + mfacb * (one - vx2) - c1o9 * oMdrho * vx2;
+ m2 = mfccb * c1o2 + mfbcb * (vvx + c1o2) + (mfacb + c1o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfacb = m0;
+ mfbcb = m1;
+ mfccb = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcac * c1o2 + mfbac * (vvx - c1o2) + (mfaac + c1o36 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcac - two* mfbac * vvx + mfaac * (one - vx2) - c1o36 * oMdrho * vx2;
+ m2 = mfcac * c1o2 + mfbac * (vvx + c1o2) + (mfaac + c1o36 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaac = m0;
+ mfbac = m1;
+ mfcac = m2;
+ ///////////c////////////////////////////////////////////////////////////////////////
+ m0 = mfcbc * c1o2 + mfbbc * (vvx - c1o2) + (mfabc + c1o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcbc - two* mfbbc * vvx + mfabc * (one - vx2) - c1o9 * oMdrho * vx2;
+ m2 = mfcbc * c1o2 + mfbbc * (vvx + c1o2) + (mfabc + c1o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfabc = m0;
+ mfbbc = m1;
+ mfcbc = m2;
+ ///////////c////////////////////////////////////////////////////////////////////////
+ m0 = mfccc * c1o2 + mfbcc * (vvx - c1o2) + (mfacc + c1o36 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfccc - two* mfbcc * vvx + mfacc * (one - vx2) - c1o36 * oMdrho * vx2;
+ m2 = mfccc * c1o2 + mfbcc * (vvx + c1o2) + (mfacc + c1o36 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfacc = m0;
+ mfbcc = m1;
+ mfccc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ (D.f[dirE])[k] = mfabb;
+ (D.f[dirW])[kw] = mfcbb;
+ (D.f[dirN])[k] = mfbab;
+ (D.f[dirS])[ks] = mfbcb;
+ (D.f[dirT])[k] = mfbba;
+ (D.f[dirB])[kb] = mfbbc;
+ (D.f[dirNE])[k] = mfaab;
+ (D.f[dirSW])[ksw] = mfccb;
+ (D.f[dirSE])[ks] = mfacb;
+ (D.f[dirNW])[kw] = mfcab;
+ (D.f[dirTE])[k] = mfaba;
+ (D.f[dirBW])[kbw] = mfcbc;
+ (D.f[dirBE])[kb] = mfabc;
+ (D.f[dirTW])[kw] = mfcba;
+ (D.f[dirTN])[k] = mfbaa;
+ (D.f[dirBS])[kbs] = mfbcc;
+ (D.f[dirBN])[kb] = mfbac;
+ (D.f[dirTS])[ks] = mfbca;
+ (D.f[dirZERO])[k] = mfbbb;
+ (D.f[dirTNE])[k] = mfaaa;
+ (D.f[dirTSE])[ks] = mfaca;
+ (D.f[dirBNE])[kb] = mfaac;
+ (D.f[dirBSE])[kbs] = mfacc;
+ (D.f[dirTNW])[kw] = mfcaa;
+ (D.f[dirTSW])[ksw] = mfcca;
+ (D.f[dirBNW])[kbw] = mfcac;
+ (D.f[dirBSW])[kbsw] = mfccc;
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ (G.g[dirE])[k] = mgabb;
+ (G.g[dirW])[kw] = mgcbb;
+ (G.g[dirN])[k] = mgbab;
+ (G.g[dirS])[ks] = mgbcb;
+ (G.g[dirT])[k] = mgbba;
+ (G.g[dirB])[kb] = mgbbc;
+ }
+ }
+}
\ No newline at end of file
diff --git a/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27_Device.cuh b/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27_Device.cuh
new file mode 100644
index 000000000..9475084ba
--- /dev/null
+++ b/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/CumulantF32018CompSP27_Device.cuh
@@ -0,0 +1,19 @@
+#ifndef LB_Kernel_Cumulant_D3Q27F3_2018_H
+#define LB_Kernel_Cumulant_D3Q27F3_2018_H
+
+#include <DataTypes.h>
+#include <curand.h>
+
+extern "C" __global__ void LB_Kernel_Cumulant_D3Q27F3_2018( real omega,
+ unsigned int* bcMatD,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ real* DDStart,
+ real* F3,
+ int size_Mat,
+ int level,
+ real* forces,
+ bool EvenOrOdd);
+
+#endif
\ No newline at end of file
diff --git a/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/package.include b/src/VirtualFluids_GPU/Kernel/Advection/Compressible/CumulantF32018/package.include
new file mode 100644
index 000000000..e69de29bb
--
GitLab