diff --git a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
index df18067a764151a1b9af50a4f2e46b3a0d46484e..de83ea55c0365b6e217cdc41a01b013fea4653da 100644
--- a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
+++ b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
@@ -161,7 +161,7 @@ void multipleLevel(const std::string& configPath)
     para->setViscosityLB(viscosityLB);
     para->setVelocityRatio( dx / dt );
     para->setViscosityRatio( dx*dx/dt );
-    para->setMainKernel("CumulantK17Almighty");
+    para->setMainKernel("CumulantK17");
 
     para->setInitialCondition([&](real coordX, real coordY, real coordZ, real &rho, real &vx, real &vy, real &vz) {
         rho = (real)0.0;
diff --git a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
index 3d6840316df78d071ca8990cc714beeba470ea23..7dc4d5096d7e3d63197a488c1c936dabb61e85ae 100644
--- a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
+++ b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
@@ -205,7 +205,7 @@ void multipleLevel(const std::string& configPath)
     bool useStreams = (nProcs > 1 ? true: false);
     // useStreams=false;
     para->setUseStreams(useStreams);
-    para->setMainKernel("CumulantK17Almighty");
+    para->setMainKernel("CumulantK17");
     para->setIsBodyForce( config.getValue<bool>("bodyForce") );
 
     para->setTimestepStartOut(uint(tStartOut/dt) );
diff --git a/src/gpu/VirtualFluids_GPU/Calculation/CollisisionStrategy.cpp b/src/gpu/VirtualFluids_GPU/Calculation/CollisisionStrategy.cpp
index 212e21f5e002201a2d4d844dd5a996a0410bc945..1c279ad523e8442bc68af079f1c789a1b130cd35 100644
--- a/src/gpu/VirtualFluids_GPU/Calculation/CollisisionStrategy.cpp
+++ b/src/gpu/VirtualFluids_GPU/Calculation/CollisisionStrategy.cpp
@@ -83,8 +83,9 @@ void CollisionAndExchange_streams::operator()(UpdateGrid27 *updateGrid, Paramete
     if (para->getUseStreams())
         para->getStreamManager()->triggerStartBulkKernel(CudaStreamIndex::Border);
 
-    //! 3. launch the collision kernel for bulk nodes
-    //!
+    //! 3. launch the collision kernel for bulk nodes. This includes nodes with \param tag Default, WriteMacroVars, ApplyBodyForce, 
+    //!    or AllFeatures. All assigned tags are listed in \param allocatedBulkFluidNodeTags during initialization in Simulation::init
+
     para->getStreamManager()->waitOnStartBulkKernelEvent(CudaStreamIndex::Bulk);
     
     for( CollisionTemplate tag: para->getParH(level)->allocatedBulkFluidNodeTags )
diff --git a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
index cb9563aafb7b019d7324ddbc31445223733b7f3a..4136614dfbfc9e0d2fc1bf7f4b01624f94eabb6f 100644
--- a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
+++ b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
@@ -55,7 +55,8 @@ void UpdateGrid27::updateGrid(int level, unsigned int t)
         refinement(this, para.get(), level);
     }
 
-
+    //////////////////////////////////////////////////////////////////////////
+    
     interactWithActuators(level, t);
 
 }
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty.cu b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17.cu
similarity index 81%
rename from src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty.cu
rename to src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17.cu
index 89b4751b8db972da94d15f8dc93dcd18cc67b2ea..54dd11142d9c16063d58330cfe7351394bdfe51c 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17.cu
@@ -1,21 +1,21 @@
-#include "CumulantK17Almighty.h"
+#include "CumulantK17.h"
 #include <logger/Logger.h>
 #include "Parameter/Parameter.h"
 #include "Parameter/CudaStreamManager.h"
-#include "CumulantK17Almighty_Device.cuh"
+#include "CumulantK17_Device.cuh"
 
 #include <cuda.h>
 
 template<TurbulenceModel turbulenceModel> 
-std::shared_ptr< CumulantK17Almighty<turbulenceModel> > CumulantK17Almighty<turbulenceModel>::getNewInstance(std::shared_ptr<Parameter> para, int level)
+std::shared_ptr< CumulantK17<turbulenceModel> > CumulantK17<turbulenceModel>::getNewInstance(std::shared_ptr<Parameter> para, int level)
 {
-	return std::shared_ptr<CumulantK17Almighty<turbulenceModel> >(new CumulantK17Almighty<turbulenceModel>(para,level));
+	return std::shared_ptr<CumulantK17<turbulenceModel> >(new CumulantK17<turbulenceModel>(para,level));
 }
 
 template<TurbulenceModel turbulenceModel>
-void CumulantK17Almighty<turbulenceModel>::run()
+void CumulantK17<turbulenceModel>::run()
 {
-	LB_Kernel_CumulantK17Almighty < turbulenceModel, false, false  > <<< cudaGrid.grid, cudaGrid.threads >>>(   para->getParD(level)->omega, 	
+	LB_Kernel_CumulantK17 < turbulenceModel, false, false  > <<< cudaGrid.grid, cudaGrid.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],	
@@ -33,18 +33,18 @@ void CumulantK17Almighty<turbulenceModel>::run()
 																												para->getParD(level)->taggedFluidNodeIndices[CollisionTemplate::Default],
         																										para->getParD(level)->numberOfTaggedFluidNodes[CollisionTemplate::Default]);
 
-	getLastCudaError("LB_Kernel_CumulantK17Almighty execution failed");
+	getLastCudaError("LB_Kernel_CumulantK17 execution failed");
 }
 
 template<TurbulenceModel turbulenceModel>
-void CumulantK17Almighty<turbulenceModel>::runOnIndices( const unsigned int *indices, unsigned int size_indices, CollisionTemplate collisionTemplate, CudaStreamIndex streamIndex )
+void CumulantK17<turbulenceModel>::runOnIndices( const unsigned int *indices, unsigned int size_indices, CollisionTemplate collisionTemplate, CudaStreamIndex streamIndex )
 {
 	cudaStream_t stream = para->getStreamManager()->getStream(streamIndex);
 	
 	switch (collisionTemplate)
 	{
 		case CollisionTemplate::Default:
-			LB_Kernel_CumulantK17Almighty < turbulenceModel, false, false  > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>(   
+			LB_Kernel_CumulantK17 < turbulenceModel, false, false  > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>(   
 																																	para->getParD(level)->omega, 	
 																																	para->getParD(level)->typeOfGridNode, 										
 																																	para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ,	
@@ -65,7 +65,7 @@ void CumulantK17Almighty<turbulenceModel>::runOnIndices( const unsigned int *ind
 			break;
 		
 		case CollisionTemplate::WriteMacroVars:
-			LB_Kernel_CumulantK17Almighty < turbulenceModel, true, false  > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>( 
+			LB_Kernel_CumulantK17 < turbulenceModel, true, false  > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>( 
 																																para->getParD(level)->omega, 	
 																																para->getParD(level)->typeOfGridNode, 										
 																																para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ,	
@@ -87,7 +87,7 @@ void CumulantK17Almighty<turbulenceModel>::runOnIndices( const unsigned int *ind
 		
 		case CollisionTemplate::Border:
 		case CollisionTemplate::AllFeatures:
-			LB_Kernel_CumulantK17Almighty < turbulenceModel, true, true  > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>(  
+			LB_Kernel_CumulantK17 < turbulenceModel, true, true  > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>(  
 																																para->getParD(level)->omega, 	
 																																para->getParD(level)->typeOfGridNode, 										
 																																para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ,	
@@ -107,7 +107,7 @@ void CumulantK17Almighty<turbulenceModel>::runOnIndices( const unsigned int *ind
 																																size_indices);
 			break;
 		case CollisionTemplate::ApplyBodyForce:
-			LB_Kernel_CumulantK17Almighty < turbulenceModel, false, true  > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>( 
+			LB_Kernel_CumulantK17 < turbulenceModel, false, true  > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>( 
 																																para->getParD(level)->omega, 	
 																																para->getParD(level)->typeOfGridNode, 										
 																																para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ,	
@@ -127,15 +127,15 @@ void CumulantK17Almighty<turbulenceModel>::runOnIndices( const unsigned int *ind
 																																size_indices);
 			break;
 		default:
-			throw std::runtime_error("Invalid CollisionTemplate in CumulantK17Almighty::runOnIndices()");
+			throw std::runtime_error("Invalid CollisionTemplate in CumulantK17::runOnIndices()");
 			break;
 	}
 
-	getLastCudaError("LB_Kernel_CumulantK17Almighty execution failed");
+	getLastCudaError("LB_Kernel_CumulantK17 execution failed");
 }
 
 template<TurbulenceModel turbulenceModel>
-CumulantK17Almighty<turbulenceModel>::CumulantK17Almighty(std::shared_ptr<Parameter> para, int level)
+CumulantK17<turbulenceModel>::CumulantK17(std::shared_ptr<Parameter> para, int level)
 {
 	this->para = para;
 	this->level = level;
@@ -150,7 +150,7 @@ CumulantK17Almighty<turbulenceModel>::CumulantK17Almighty(std::shared_ptr<Parame
 	VF_LOG_INFO("Using turbulence model: {}", turbulenceModel);
 }
 
-template class CumulantK17Almighty<TurbulenceModel::AMD>;
-template class CumulantK17Almighty<TurbulenceModel::Smagorinsky>;
-template class CumulantK17Almighty<TurbulenceModel::QR>;
-template class CumulantK17Almighty<TurbulenceModel::None>;
+template class CumulantK17<TurbulenceModel::AMD>;
+template class CumulantK17<TurbulenceModel::Smagorinsky>;
+template class CumulantK17<TurbulenceModel::QR>;
+template class CumulantK17<TurbulenceModel::None>;
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17.h b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17.h
new file mode 100644
index 0000000000000000000000000000000000000000..00c79a30c9ccf9a89901165d020fc85d5a479c1d
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17.h
@@ -0,0 +1,20 @@
+#ifndef CUMULANT_K17_H
+#define CUMULANT_K17_H
+
+#include "Kernel/KernelImp.h"
+#include "Parameter/Parameter.h"
+
+template<TurbulenceModel turbulenceModel> 
+class CumulantK17 : public KernelImp
+{
+public:
+	static std::shared_ptr< CumulantK17<turbulenceModel> > getNewInstance(std::shared_ptr< Parameter> para, int level);
+	void run() override;
+    void runOnIndices(const unsigned int *indices, unsigned int size_indices, CollisionTemplate collisionTemplate, CudaStreamIndex streamIndex) override;
+
+private:
+    CumulantK17();
+    CumulantK17(std::shared_ptr<Parameter> para, int level);
+};
+
+#endif 
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp.cu b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp.cu
deleted file mode 100644
index b176b94d07e7f280d738a797d5bd853095e3caed..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp.cu
+++ /dev/null
@@ -1,33 +0,0 @@
-#include "CumulantK17Comp.h"
-
-#include "Parameter/Parameter.h"
-#include "CumulantK17Comp_Device.cuh"
-#include "cuda/CudaGrid.h"
-
-std::shared_ptr<CumulantK17Comp> CumulantK17Comp::getNewInstance(std::shared_ptr<Parameter> para, int level)
-{
-	return std::shared_ptr<CumulantK17Comp>(new CumulantK17Comp(para,level));
-}
-
-void CumulantK17Comp::run()
-{
-	LB_Kernel_CumulantK17Comp <<< cudaGrid.grid, cudaGrid.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)->numberOfNodes,
-													level,
-													para->getForcesDev(),
-                                                    para->getQuadricLimitersDev(),
-													para->getParD(level)->isEvenTimestep);
-	getLastCudaError("LB_Kernel_CumulantK17Comp execution failed");
-}
-
-CumulantK17Comp::CumulantK17Comp(std::shared_ptr<Parameter> para, int level): KernelImp(para, level)
-{
-	myPreProcessorTypes.push_back(InitCompSP27);
-	myKernelGroup = BasicKernel;
-	this->cudaGrid = vf::cuda::CudaGrid(para->getParD(level)->numberofthreads, para->getParD(level)->numberOfNodes);
-}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp.h b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp.h
deleted file mode 100644
index 22a95a688e5d078d7b710f494bfea360c9af0d6b..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp.h
+++ /dev/null
@@ -1,17 +0,0 @@
-#ifndef CUMULANT_K17_COMP_H
-#define CUMULANT_K17_COMP_H
-
-#include "Kernel/KernelImp.h"
-
-class CumulantK17Comp : public KernelImp
-{
-public:
-	static std::shared_ptr<CumulantK17Comp> getNewInstance(std::shared_ptr< Parameter> para, int level);
-	void run();
-
-private:
-	CumulantK17Comp();
-	CumulantK17Comp(std::shared_ptr< Parameter> para, int level);
-};
-
-#endif 
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp_Device.cu b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp_Device.cu
deleted file mode 100644
index 7cf27aa883cbfd3a0e4a0a36fa61649a62d06eeb..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp_Device.cu
+++ /dev/null
@@ -1,1040 +0,0 @@
-#include "LBM/LB.h" 
-#include "lbm/constants/D3Q27.h"
-#include <lbm/constants/NumericConstants.h>
-
-using namespace vf::lbm::constant;
-using namespace vf::lbm::dir;
-#include "math.h"
-
-
-__global__ void LB_Kernel_CumulantK17Comp(real omega,
-	unsigned int* bcMatD,
-	unsigned int* neighborX,
-	unsigned int* neighborY,
-	unsigned int* neighborZ,
-	real* DDStart,
-	int size_Mat,
-	int level,
-	real* forces,
-	real* quadricLimiters,
-	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[DIR_P00] = &DDStart[DIR_P00   *size_Mat];
-				D.f[DIR_M00] = &DDStart[DIR_M00   *size_Mat];
-				D.f[DIR_0P0] = &DDStart[DIR_0P0   *size_Mat];
-				D.f[DIR_0M0] = &DDStart[DIR_0M0   *size_Mat];
-				D.f[DIR_00P] = &DDStart[DIR_00P   *size_Mat];
-				D.f[DIR_00M] = &DDStart[DIR_00M   *size_Mat];
-				D.f[DIR_PP0] = &DDStart[DIR_PP0  *size_Mat];
-				D.f[DIR_MM0] = &DDStart[DIR_MM0  *size_Mat];
-				D.f[DIR_PM0] = &DDStart[DIR_PM0  *size_Mat];
-				D.f[DIR_MP0] = &DDStart[DIR_MP0  *size_Mat];
-				D.f[DIR_P0P] = &DDStart[DIR_P0P  *size_Mat];
-				D.f[DIR_M0M] = &DDStart[DIR_M0M  *size_Mat];
-				D.f[DIR_P0M] = &DDStart[DIR_P0M  *size_Mat];
-				D.f[DIR_M0P] = &DDStart[DIR_M0P  *size_Mat];
-				D.f[DIR_0PP] = &DDStart[DIR_0PP  *size_Mat];
-				D.f[DIR_0MM] = &DDStart[DIR_0MM  *size_Mat];
-				D.f[DIR_0PM] = &DDStart[DIR_0PM  *size_Mat];
-				D.f[DIR_0MP] = &DDStart[DIR_0MP  *size_Mat];
-				D.f[DIR_000] = &DDStart[DIR_000*size_Mat];
-				D.f[DIR_PPP] = &DDStart[DIR_PPP *size_Mat];
-				D.f[DIR_MMP] = &DDStart[DIR_MMP *size_Mat];
-				D.f[DIR_PMP] = &DDStart[DIR_PMP *size_Mat];
-				D.f[DIR_MPP] = &DDStart[DIR_MPP *size_Mat];
-				D.f[DIR_PPM] = &DDStart[DIR_PPM *size_Mat];
-				D.f[DIR_MMM] = &DDStart[DIR_MMM *size_Mat];
-				D.f[DIR_PMM] = &DDStart[DIR_PMM *size_Mat];
-				D.f[DIR_MPM] = &DDStart[DIR_MPM *size_Mat];
-			}
-			else
-			{
-				D.f[DIR_M00] = &DDStart[DIR_P00   *size_Mat];
-				D.f[DIR_P00] = &DDStart[DIR_M00   *size_Mat];
-				D.f[DIR_0M0] = &DDStart[DIR_0P0   *size_Mat];
-				D.f[DIR_0P0] = &DDStart[DIR_0M0   *size_Mat];
-				D.f[DIR_00M] = &DDStart[DIR_00P   *size_Mat];
-				D.f[DIR_00P] = &DDStart[DIR_00M   *size_Mat];
-				D.f[DIR_MM0] = &DDStart[DIR_PP0  *size_Mat];
-				D.f[DIR_PP0] = &DDStart[DIR_MM0  *size_Mat];
-				D.f[DIR_MP0] = &DDStart[DIR_PM0  *size_Mat];
-				D.f[DIR_PM0] = &DDStart[DIR_MP0  *size_Mat];
-				D.f[DIR_M0M] = &DDStart[DIR_P0P  *size_Mat];
-				D.f[DIR_P0P] = &DDStart[DIR_M0M  *size_Mat];
-				D.f[DIR_M0P] = &DDStart[DIR_P0M  *size_Mat];
-				D.f[DIR_P0M] = &DDStart[DIR_M0P  *size_Mat];
-				D.f[DIR_0MM] = &DDStart[DIR_0PP  *size_Mat];
-				D.f[DIR_0PP] = &DDStart[DIR_0MM  *size_Mat];
-				D.f[DIR_0MP] = &DDStart[DIR_0PM  *size_Mat];
-				D.f[DIR_0PM] = &DDStart[DIR_0MP  *size_Mat];
-				D.f[DIR_000] = &DDStart[DIR_000*size_Mat];
-				D.f[DIR_MMM] = &DDStart[DIR_PPP *size_Mat];
-				D.f[DIR_PPM] = &DDStart[DIR_MMP *size_Mat];
-				D.f[DIR_MPM] = &DDStart[DIR_PMP *size_Mat];
-				D.f[DIR_PMM] = &DDStart[DIR_MPP *size_Mat];
-				D.f[DIR_MMP] = &DDStart[DIR_PPM *size_Mat];
-				D.f[DIR_PPP] = &DDStart[DIR_MMM *size_Mat];
-				D.f[DIR_MPP] = &DDStart[DIR_PMM *size_Mat];
-				D.f[DIR_PMP] = &DDStart[DIR_MPM *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];
-
-			//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 mfcbb = (D.f[DIR_P00])[k];//[ke   ];// +  c2over27 ;(D.f[DIR_P00   ])[k  ];//ke
-			real mfabb = (D.f[DIR_M00])[kw];//[kw   ];// +  c2over27 ;(D.f[DIR_M00   ])[kw ];
-			real mfbcb = (D.f[DIR_0P0])[k];//[kn   ];// +  c2over27 ;(D.f[DIR_0P0   ])[k  ];//kn
-			real mfbab = (D.f[DIR_0M0])[ks];//[ks   ];// +  c2over27 ;(D.f[DIR_0M0   ])[ks ];
-			real mfbbc = (D.f[DIR_00P])[k];//[kt   ];// +  c2over27 ;(D.f[DIR_00P   ])[k  ];//kt
-			real mfbba = (D.f[DIR_00M])[kb];//[kb   ];// +  c2over27 ;(D.f[DIR_00M   ])[kb ];
-			real mfccb = (D.f[DIR_PP0])[k];//[kne  ];// +  c1over54 ;(D.f[DIR_PP0  ])[k  ];//kne
-			real mfaab = (D.f[DIR_MM0])[ksw];//[ksw  ];// +  c1over54 ;(D.f[DIR_MM0  ])[ksw];
-			real mfcab = (D.f[DIR_PM0])[ks];//[kse  ];// +  c1over54 ;(D.f[DIR_PM0  ])[ks ];//kse
-			real mfacb = (D.f[DIR_MP0])[kw];//[knw  ];// +  c1over54 ;(D.f[DIR_MP0  ])[kw ];//knw
-			real mfcbc = (D.f[DIR_P0P])[k];//[kte  ];// +  c1over54 ;(D.f[DIR_P0P  ])[k  ];//kte
-			real mfaba = (D.f[DIR_M0M])[kbw];//[kbw  ];// +  c1over54 ;(D.f[DIR_M0M  ])[kbw];
-			real mfcba = (D.f[DIR_P0M])[kb];//[kbe  ];// +  c1over54 ;(D.f[DIR_P0M  ])[kb ];//kbe
-			real mfabc = (D.f[DIR_M0P])[kw];//[ktw  ];// +  c1over54 ;(D.f[DIR_M0P  ])[kw ];//ktw
-			real mfbcc = (D.f[DIR_0PP])[k];//[ktn  ];// +  c1over54 ;(D.f[DIR_0PP  ])[k  ];//ktn
-			real mfbaa = (D.f[DIR_0MM])[kbs];//[kbs  ];// +  c1over54 ;(D.f[DIR_0MM  ])[kbs];
-			real mfbca = (D.f[DIR_0PM])[kb];//[kbn  ];// +  c1over54 ;(D.f[DIR_0PM  ])[kb ];//kbn
-			real mfbac = (D.f[DIR_0MP])[ks];//[kts  ];// +  c1over54 ;(D.f[DIR_0MP  ])[ks ];//kts
-			real mfbbb = (D.f[DIR_000])[k];//[kzero];// +  c8over27 ;(D.f[DIR_000])[k  ];//kzero
-			real mfccc = (D.f[DIR_PPP])[k];//[ktne ];// +  c1over216;(D.f[DIR_PPP ])[k  ];//ktne
-			real mfaac = (D.f[DIR_MMP])[ksw];//[ktsw ];// +  c1over216;(D.f[DIR_MMP ])[ksw];//ktsw
-			real mfcac = (D.f[DIR_PMP])[ks];//[ktse ];// +  c1over216;(D.f[DIR_PMP ])[ks ];//ktse
-			real mfacc = (D.f[DIR_MPP])[kw];//[ktnw ];// +  c1over216;(D.f[DIR_MPP ])[kw ];//ktnw
-			real mfcca = (D.f[DIR_PPM])[kb];//[kbne ];// +  c1over216;(D.f[DIR_PPM ])[kb ];//kbne
-			real mfaaa = (D.f[DIR_MMM])[kbsw];//[kbsw ];// +  c1over216;(D.f[DIR_MMM ])[kbsw];
-			real mfcaa = (D.f[DIR_PMM])[kbs];//[kbse ];// +  c1over216;(D.f[DIR_PMM ])[kbs];//kbse
-			real mfaca = (D.f[DIR_MPM])[kbw];//[kbnw ];// +  c1over216;(D.f[DIR_MPM ])[kbw];//kbnw
-											////////////////////////////////////////////////////////////////////////////////////
-			real drho = ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) +
-				(((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
-				((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb;
-
-			real rho = c1o1 + drho;
-			////////////////////////////////////////////////////////////////////////////////////
-			//slow
-			//real oMdrho = one - ((((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)));//fehlt mfbbb
-			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((double)c2o1, (double)level)); //zero;//0.0032653/(pow(two,level)); //0.000000005;//(two/1600000.0) / 120.0; //
-			real fy = forces[1] / (pow((double)c2o1, (double)level)); //zero;
-			real fz = forces[2] / (pow((double)c2o1, (double)level)); //zero;
-			vvx += fx*c1o2;
-			vvy += fy*c1o2;
-			vvz += fz*c1o2;
-			////////////////////////////////////////////////////////////////////////////////////
-			//real omega = omega_in;
-			////////////////////////////////////////////////////////////////////////////////////
-			//fast
-			real oMdrho = c1o1; // comp special
-							   //real oMdrho = one - (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);//fehlt mfbbb nicht mehr
-							   //real vvx    =mfccc-mfaaa + mfcac-mfaca + mfcaa-mfacc + mfcca-mfaac + 
-							   //				mfcba-mfabc + mfcbc-mfaba + mfcab-mfacb + mfccb-mfaab +
-							   //				mfcbb-mfabb;
-							   //real vvy    =mfccc-mfaaa + mfaca-mfcac + mfacc-mfcaa + mfcca-mfaac + 
-							   //				mfbca-mfbac + mfbcc-mfbaa + mfacb-mfcab + mfccb-mfaab +
-							   //				mfbcb-mfbab;
-							   //real vvz    =mfccc-mfaaa + mfcac-mfaca + mfacc-mfcaa + mfaac-mfcca + 
-							   //				mfbac-mfbca + mfbcc-mfbaa + mfabc-mfcba + mfcbc-mfaba +
-							   //				mfbbc-mfbba;
-							   ////////////////////////////////////////////////////////////////////////////////////
-							   // oMdrho assembler style -------> faaaaaastaaaa
-							   // or much sloooowaaaa ... it dep�ndssssss on sadaku
-			real m0, m1, m2;
-			//real oMdrho;
-			//{
-			//	oMdrho=mfccc+mfaaa;
-			//	m0=mfaca+mfcac;
-			//	m1=mfacc+mfcaa;
-			//	m2=mfaac+mfcca;
-			//	oMdrho+=m0;
-			//	m1+=m2;
-			//	oMdrho+=m1;
-			//	m0=mfbac+mfbca;
-			//	m1=mfbaa+mfbcc;
-			//	m0+=m1;
-			//	m1=mfabc+mfcba;
-			//	m2=mfaba+mfcbc;
-			//	m1+=m2;
-			//	m0+=m1;
-			//	m1=mfacb+mfcab;
-			//	m2=mfaab+mfccb;
-			//	m1+=m2;
-			//	m0+=m1;
-			//	oMdrho+=m0;
-			//	m0=mfabb+mfcbb;
-			//	m1=mfbab+mfbcb;
-			//	m2=mfbba+mfbbc;
-			//	m0+=m1+m2;
-			//	m0+=mfbbb; //hat gefehlt
-			//	oMdrho = one - (oMdrho + m0);
-			//}
-			//real vvx;
-			real vx2;
-			//{
-			//	vvx = mfccc-mfaaa;
-			//	m0  = mfcac-mfaca;
-			//	m1  = mfcaa-mfacc;
-			//	m2  = mfcca-mfaac;
-			//	vvx+= m0;
-			//	m1 += m2;
-			//	vvx+= m1;
-			//	vx2 = mfcba-mfabc;
-			//	m0  = mfcbc-mfaba;
-			//	m1  = mfcab-mfacb;
-			//	m2  = mfccb-mfaab;
-			//	vx2+= m0;
-			//	m1 += m2;
-			//	vx2+= m1;
-			//	vvx+= vx2;
-			//	vx2 = mfcbb-mfabb;
-			//	vvx+= vx2;
-			//}
-			//real vvy;
-			real vy2;
-			//{
-			//	vvy = mfccc-mfaaa;
-			//	m0  = mfaca-mfcac;
-			//	m1  = mfacc-mfcaa;
-			//	m2  = mfcca-mfaac;
-			//	vvy+= m0;
-			//	m1 += m2;
-			//	vvy+= m1;
-			//	vy2 = mfbca-mfbac;
-			//	m0  = mfbcc-mfbaa;
-			//	m1  = mfacb-mfcab;
-			//	m2  = mfccb-mfaab;
-			//	vy2+= m0;
-			//	m1 += m2;
-			//	vy2+= m1;
-			//	vvy+= vy2;
-			//	vy2 = mfbcb-mfbab;
-			//	vvy+= vy2;
-			//}
-			//real vvz;
-			real vz2;
-			//{
-			//	vvz = mfccc-mfaaa;
-			//	m0  = mfcac-mfaca;
-			//	m1  = mfacc-mfcaa;
-			//	m2  = mfaac-mfcca;
-			//	vvz+= m0;
-			//	m1 += m2;
-			//	vvz+= m1;
-			//	vz2 = mfbac-mfbca;
-			//	m0  = mfbcc-mfbaa;
-			//	m1  = mfabc-mfcba;
-			//	m2  = mfcbc-mfaba;
-			//	vz2+= m0;
-			//	m1 += m2;
-			//	vz2+= m1;
-			//	vvz+= vz2;
-			//	vz2 = mfbbc-mfbba;
-			//	vvz+= vz2;
-			//}
-			vx2 = vvx*vvx;
-			vy2 = vvy*vvy;
-			vz2 = vvz*vvz;
-			////////////////////////////////////////////////////////////////////////////////////
-			real wadjust;
-			real qudricLimitP = quadricLimiters[0];  //0.01f; //  * 0.0001f; // 1000000.0f; // 1000000.0f; //
-			real qudricLimitM = quadricLimiters[1];  //0.01f; //  * 0.0001f; // 1000000.0f; // 1000000.0f; //
-			real qudricLimitD = quadricLimiters[2];  //0.01f; //  * 0.001f;  // 1000000.0f; // 1000000.0f; //
-									  ////////////////////////////////////////////////////////////////////////////////////
-									  //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 - c2o1*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfaba + mfabc;
-			m1 = mfabc - mfaba;
-			m0 = m2 + mfabb;
-			mfaba = m0;
-			m0 += c1o9 * oMdrho;
-			mfabb = m1 - m0 * vvz;
-			mfabc = m2 - c2o1*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfaca + mfacc;
-			m1 = mfacc - mfaca;
-			m0 = m2 + mfacb;
-			mfaca = m0;
-			m0 += c1o36 * oMdrho;
-			mfacb = m1 - m0 * vvz;
-			mfacc = m2 - c2o1*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfbaa + mfbac;
-			m1 = mfbac - mfbaa;
-			m0 = m2 + mfbab;
-			mfbaa = m0;
-			m0 += c1o9 * oMdrho;
-			mfbab = m1 - m0 * vvz;
-			mfbac = m2 - c2o1*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfbba + mfbbc;
-			m1 = mfbbc - mfbba;
-			m0 = m2 + mfbbb;
-			mfbba = m0;
-			m0 += c4o9 * oMdrho;
-			mfbbb = m1 - m0 * vvz;
-			mfbbc = m2 - c2o1*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfbca + mfbcc;
-			m1 = mfbcc - mfbca;
-			m0 = m2 + mfbcb;
-			mfbca = m0;
-			m0 += c1o9 * oMdrho;
-			mfbcb = m1 - m0 * vvz;
-			mfbcc = m2 - c2o1*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfcaa + mfcac;
-			m1 = mfcac - mfcaa;
-			m0 = m2 + mfcab;
-			mfcaa = m0;
-			m0 += c1o36 * oMdrho;
-			mfcab = m1 - m0 * vvz;
-			mfcac = m2 - c2o1*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfcba + mfcbc;
-			m1 = mfcbc - mfcba;
-			m0 = m2 + mfcbb;
-			mfcba = m0;
-			m0 += c1o9 * oMdrho;
-			mfcbb = m1 - m0 * vvz;
-			mfcbc = m2 - c2o1*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfcca + mfccc;
-			m1 = mfccc - mfcca;
-			m0 = m2 + mfccb;
-			mfcca = m0;
-			m0 += c1o36 * oMdrho;
-			mfccb = m1 - m0 * vvz;
-			mfccc = m2 - c2o1*	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 - c2o1*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfaab + mfacb;
-			m1 = mfacb - mfaab;
-			m0 = m2 + mfabb;
-			mfaab = m0;
-			mfabb = m1 - m0 * vvy;
-			mfacb = m2 - c2o1*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfaac + mfacc;
-			m1 = mfacc - mfaac;
-			m0 = m2 + mfabc;
-			mfaac = m0;
-			m0 += c1o18 * oMdrho;
-			mfabc = m1 - m0 * vvy;
-			mfacc = m2 - c2o1*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfbaa + mfbca;
-			m1 = mfbca - mfbaa;
-			m0 = m2 + mfbba;
-			mfbaa = m0;
-			m0 += c2o3 * oMdrho;
-			mfbba = m1 - m0 * vvy;
-			mfbca = m2 - c2o1*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfbab + mfbcb;
-			m1 = mfbcb - mfbab;
-			m0 = m2 + mfbbb;
-			mfbab = m0;
-			mfbbb = m1 - m0 * vvy;
-			mfbcb = m2 - c2o1*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfbac + mfbcc;
-			m1 = mfbcc - mfbac;
-			m0 = m2 + mfbbc;
-			mfbac = m0;
-			m0 += c2o9 * oMdrho;
-			mfbbc = m1 - m0 * vvy;
-			mfbcc = m2 - c2o1*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfcaa + mfcca;
-			m1 = mfcca - mfcaa;
-			m0 = m2 + mfcba;
-			mfcaa = m0;
-			m0 += c1o6 * oMdrho;
-			mfcba = m1 - m0 * vvy;
-			mfcca = m2 - c2o1*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfcab + mfccb;
-			m1 = mfccb - mfcab;
-			m0 = m2 + mfcbb;
-			mfcab = m0;
-			mfcbb = m1 - m0 * vvy;
-			mfccb = m2 - c2o1*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfcac + mfccc;
-			m1 = mfccc - mfcac;
-			m0 = m2 + mfcbc;
-			mfcac = m0;
-			m0 += c1o18 * oMdrho;
-			mfcbc = m1 - m0 * vvy;
-			mfccc = m2 - c2o1*	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 += c1o1* oMdrho;
-			mfbaa = m1 - m0 * vvx;
-			mfcaa = m2 - c2o1*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfaba + mfcba;
-			m1 = mfcba - mfaba;
-			m0 = m2 + mfbba;
-			mfaba = m0;
-			mfbba = m1 - m0 * vvx;
-			mfcba = m2 - c2o1*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfaca + mfcca;
-			m1 = mfcca - mfaca;
-			m0 = m2 + mfbca;
-			mfaca = m0;
-			m0 += c1o3 * oMdrho;
-			mfbca = m1 - m0 * vvx;
-			mfcca = m2 - c2o1*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfaab + mfcab;
-			m1 = mfcab - mfaab;
-			m0 = m2 + mfbab;
-			mfaab = m0;
-			mfbab = m1 - m0 * vvx;
-			mfcab = m2 - c2o1*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfabb + mfcbb;
-			m1 = mfcbb - mfabb;
-			m0 = m2 + mfbbb;
-			mfabb = m0;
-			mfbbb = m1 - m0 * vvx;
-			mfcbb = m2 - c2o1*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfacb + mfccb;
-			m1 = mfccb - mfacb;
-			m0 = m2 + mfbcb;
-			mfacb = m0;
-			mfbcb = m1 - m0 * vvx;
-			mfccb = m2 - c2o1*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfaac + mfcac;
-			m1 = mfcac - mfaac;
-			m0 = m2 + mfbac;
-			mfaac = m0;
-			m0 += c1o3 * oMdrho;
-			mfbac = m1 - m0 * vvx;
-			mfcac = m2 - c2o1*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfabc + mfcbc;
-			m1 = mfcbc - mfabc;
-			m0 = m2 + mfbbc;
-			mfabc = m0;
-			mfbbc = m1 - m0 * vvx;
-			mfcbc = m2 - c2o1*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2 = mfacc + mfccc;
-			m1 = mfccc - mfacc;
-			m0 = m2 + mfbcc;
-			mfacc = m0;
-			m0 += c1o9 * oMdrho;
-			mfbcc = m1 - m0 * vvx;
-			mfccc = m2 - c2o1*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-
-			////////////////////////////////////////////////////////////////////////////////////
-			// Cumulants
-			////////////////////////////////////////////////////////////////////////////////////
-			real OxxPyyPzz = c1o1;	//set the bulk viscosity one is high / two is very low and zero is (too) high ... (also called omega 2)
-
-									////////////////////////////////////////////////////////////
-									//3.
-									//////////////////////////////
-			real OxyyPxzz = c8o1*(-c2o1 + omega)*(c1o1 + c2o1*omega) / (-c8o1 - c14o1*omega + c7o1*omega*omega);//one;
-			real OxyyMxzz = c8o1*(-c2o1 + omega)*(-c7o1 + c4o1*omega) / (c56o1 - c50o1*omega + c9o1*omega*omega);//one;
-			real Oxyz = c24o1*(-c2o1 + omega)*(-c2o1 - c7o1*omega + c3o1*omega*omega) / (c48o1 + c152o1*omega - c130o1*omega*omega + c29o1*omega*omega*omega);//one;
-																																										  ////////////////////////////////////////////////////////////
-																																										  //4.
-																																										  //////////////////////////////
-			real O4 = c1o1;
-			//////////////////////////////
-			//real O4        = omega;//TRT
-			////////////////////////////////////////////////////////////
-			//5.
-			//////////////////////////////
-			real O5 = c1o1;
-			////////////////////////////////////////////////////////////
-			//6.
-			//////////////////////////////
-			real O6 = c1o1;
-			////////////////////////////////////////////////////////////
-
-
-			//central moments to cumulants
-			//4.
-			real CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + c2o1 * mfbba * mfbab) / rho;
-			real CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + c2o1 * mfbba * mfabb) / rho;
-			real CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + c2o1 * mfbab * mfabb) / rho;
-
-			real CUMcca = mfcca - (((mfcaa * mfaca + c2o1 * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho - c1o9*(drho / rho));
-			real CUMcac = mfcac - (((mfcaa * mfaac + c2o1 * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho - c1o9*(drho / rho));
-			real CUMacc = mfacc - (((mfaac * mfaca + c2o1 * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho - c1o9*(drho / rho));
-
-			//5.
-			real CUMbcc = mfbcc - ((mfaac * mfbca + mfaca * mfbac + c4o1 * mfabb * mfbbb + c2o1 * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac)) / rho;
-			real CUMcbc = mfcbc - ((mfaac * mfcba + mfcaa * mfabc + c4o1 * mfbab * mfbbb + c2o1 * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc)) / rho;
-			real CUMccb = mfccb - ((mfcaa * mfacb + mfaca * mfcab + c4o1 * mfbba * mfbbb + c2o1 * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab)) / rho;
-
-			//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)) / rho
-				+ (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
-					+ c2o1 * (mfcaa * mfaca * mfaac)
-					+ c16o1 *  mfbba * mfbab * mfabb) / (rho * rho)
-				- c1o3 * (mfacc + mfcac + mfcca) / rho
-				- c1o9 * (mfcaa + mfaca + mfaac) / rho
-				+ (c2o1 * (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 = -c3o1*omega*mfbba;
-			real Dxz = -c3o1*omega*mfbab;
-			real Dyz = -c3o1*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;
-
-			//relax
-			mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz) - c3o1 * (c1o1 - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);//-magicBulk*OxxPyyPzz;
-			mxxMyy += omega * (-mxxMyy) - c3o1 * (c1o1 + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
-			mxxMzz += omega * (-mxxMzz) - c3o1 * (c1o1 + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
-
-			///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-			/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-			////no correction
-			//mxxPyyPzz += OxxPyyPzz*(mfaaa-mxxPyyPzz);//-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 * (-c2o1*  mxxMyy + mxxMzz + mxxPyyPzz);
-			mfaac = c1o3 * (mxxMyy - c2o1* mxxMzz + mxxPyyPzz);
-
-
-			//relax
-			//////////////////////////////////////////////////////////////////////////
-			//das ist der limiter
-			wadjust = Oxyz + (c1o1 - Oxyz)*abs(mfbbb) / (abs(mfbbb) + qudricLimitD);
-			mfbbb += wadjust * (-mfbbb);
-			wadjust = OxyyPxzz + (c1o1 - OxyyPxzz)*abs(mxxyPyzz) / (abs(mxxyPyzz) + qudricLimitP);
-			mxxyPyzz += wadjust * (-mxxyPyzz);
-			wadjust = OxyyMxzz + (c1o1 - OxyyMxzz)*abs(mxxyMyzz) / (abs(mxxyMyzz) + qudricLimitM);
-			mxxyMyzz += wadjust * (-mxxyMyzz);
-			wadjust = OxyyPxzz + (c1o1 - OxyyPxzz)*abs(mxxzPyyz) / (abs(mxxzPyyz) + qudricLimitP);
-			mxxzPyyz += wadjust * (-mxxzPyyz);
-			wadjust = OxyyMxzz + (c1o1 - OxyyMxzz)*abs(mxxzMyyz) / (abs(mxxzMyyz) + qudricLimitM);
-			mxxzMyyz += wadjust * (-mxxzMyyz);
-			wadjust = OxyyPxzz + (c1o1 - OxyyPxzz)*abs(mxyyPxzz) / (abs(mxyyPxzz) + qudricLimitP);
-			mxyyPxzz += wadjust * (-mxyyPxzz);
-			wadjust = OxyyMxzz + (c1o1 - OxyyMxzz)*abs(mxyyMxzz) / (abs(mxyyMxzz) + qudricLimitM);
-			mxyyMxzz += wadjust * (-mxyyMxzz);
-			//////////////////////////////////////////////////////////////////////////
-			//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); 
-			//////////////////////////////////////////////////////////////////////////
-			real factorA = (c4o1 + c2o1*omega - c3o1*omega*omega) / (c2o1 - c7o1*omega + c5o1*omega*omega);
-			real factorB = (c4o1 + c28o1*omega - c14o1*omega*omega) / (c6o1 - c21o1*omega + c15o1*omega*omega);
-			//////////////////////////////////////////////////////////////////////////
-			//ohne limiter
-			//CUMacc += O4 * (-CUMacc); 
-			//CUMcac += O4 * (-CUMcac); 
-			//CUMcca += O4 * (-CUMcca); 
-			//CUMbbc += O4 * (-CUMbbc); 
-			//CUMbcb += O4 * (-CUMbcb); 
-			//CUMcbb += O4 * (-CUMcbb); 
-			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);
-			//////////////////////////////////////////////////////////////////////////
-
-
-			//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 + c2o1 * mfbba * mfbab) / rho;
-			mfbcb = CUMbcb + ((mfaca + c1o3) * mfbab + c2o1 * mfbba * mfabb) / rho;
-			mfbbc = CUMbbc + ((mfaac + c1o3) * mfbba + c2o1 * mfbab * mfabb) / rho;
-
-			mfcca = CUMcca + (((mfcaa * mfaca + c2o1 * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho - c1o9*(drho / rho));
-			mfcac = CUMcac + (((mfcaa * mfaac + c2o1 * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho - c1o9*(drho / rho));
-			mfacc = CUMacc + (((mfaac * mfaca + c2o1 * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho - c1o9*(drho / rho));
-
-			//5.
-			mfbcc = CUMbcc + ((mfaac * mfbca + mfaca * mfbac + c4o1 * mfabb * mfbbb + c2o1 * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac)) / rho;
-			mfcbc = CUMcbc + ((mfaac * mfcba + mfcaa * mfabc + c4o1 * mfbab * mfbbb + c2o1 * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc)) / rho;
-			mfccb = CUMccb + ((mfcaa * mfacb + mfaca * mfcab + c4o1 * mfbba * mfbbb + c2o1 * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab)) / rho;
-
-			//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)) / rho
-				+ (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
-					+ c2o1 * (mfcaa * mfaca * mfaac)
-					+ c16o1 *  mfbba * mfbab * mfabb) / (rho * rho)
-				- c1o3 * (mfacc + mfcac + mfcca) / rho
-				- c1o9 * (mfcaa + mfaca + mfaac) / rho
-				+ (c2o1 * (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 + c1o1* oMdrho) * (vz2 - vvz) * c1o2;
-			m1 = -mfaac - c2o1* mfaab *  vvz + mfaaa                * (c1o1 - vz2) - c1o1* oMdrho * vz2;
-			m2 = mfaac * c1o2 + mfaab * (vvz + c1o2) + (mfaaa + c1o1* oMdrho) * (vz2 + vvz) * c1o2;
-			mfaaa = m0;
-			mfaab = m1;
-			mfaac = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 = mfabc * c1o2 + mfabb * (vvz - c1o2) + mfaba * (vz2 - vvz) * c1o2;
-			m1 = -mfabc - c2o1* mfabb *  vvz + mfaba * (c1o1 - 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 - c2o1* mfacb *  vvz + mfaca                  * (c1o1 - 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 - c2o1* mfbab *  vvz + mfbaa * (c1o1 - 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 - c2o1* mfbbb *  vvz + mfbba * (c1o1 - 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 - c2o1* mfbcb *  vvz + mfbca * (c1o1 - 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 - c2o1* mfcab *  vvz + mfcaa                  * (c1o1 - 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 - c2o1* mfcbb *  vvz + mfcba * (c1o1 - 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 - c2o1* mfccb *  vvz + mfcca                  * (c1o1 - 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 - c2o1* mfaba *  vvy + mfaaa                  * (c1o1 - 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 - c2o1* mfabb *  vvy + mfaab                  * (c1o1 - 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 - c2o1* mfabc *  vvy + mfaac                  * (c1o1 - 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 - c2o1* mfbba *  vvy + mfbaa * (c1o1 - 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 - c2o1* mfbbb *  vvy + mfbab * (c1o1 - 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 - c2o1* mfbbc *  vvy + mfbac * (c1o1 - 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 - c2o1* mfcba *  vvy + mfcaa                   * (c1o1 - 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 - c2o1* mfcbb *  vvy + mfcab                  * (c1o1 - 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 - c2o1* mfcbc *  vvy + mfcac                   * (c1o1 - 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 - c2o1* mfbaa *  vvx + mfaaa                   * (c1o1 - 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 - c2o1* mfbba *  vvx + mfaba                  * (c1o1 - 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 - c2o1* mfbca *  vvx + mfaca                   * (c1o1 - 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 - c2o1* mfbab *  vvx + mfaab                  * (c1o1 - 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 - c2o1* mfbbb *  vvx + mfabb                  * (c1o1 - 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 - c2o1* mfbcb *  vvx + mfacb                  * (c1o1 - 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 - c2o1* mfbac *  vvx + mfaac                   * (c1o1 - 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 - c2o1* mfbbc *  vvx + mfabc                  * (c1o1 - 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 - c2o1* mfbcc *  vvx + mfacc                   * (c1o1 - vx2) - c1o36 * oMdrho * vx2;
-			m2 = mfccc * c1o2 + mfbcc * (vvx + c1o2) + (mfacc + c1o36 * oMdrho) * (vx2 + vvx) * c1o2;
-			mfacc = m0;
-			mfbcc = m1;
-			mfccc = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-
-			////////////////////////////////////////////////////////////////////////////////////
-			(D.f[DIR_P00])[k] = mfabb;//(D.f[ DIR_P00   ])[ke   ] = mfabb;// -  c2over27 ;  (D.f[ DIR_P00   ])[k   ]                                                                     
-			(D.f[DIR_M00])[kw] = mfcbb;//(D.f[ DIR_M00   ])[kw   ] = mfcbb;// -  c2over27 ;  (D.f[ DIR_M00   ])[kw  ]                                                                   
-			(D.f[DIR_0P0])[k] = mfbab;//(D.f[ DIR_0P0   ])[kn   ] = mfbab;// -  c2over27 ;	 (D.f[ DIR_0P0   ])[k   ]
-			(D.f[DIR_0M0])[ks] = mfbcb;//(D.f[ DIR_0M0   ])[ks   ] = mfbcb;// -  c2over27 ;	 (D.f[ DIR_0M0   ])[ks  ]
-			(D.f[DIR_00P])[k] = mfbba;//(D.f[ DIR_00P   ])[kt   ] = mfbba;// -  c2over27 ;	 (D.f[ DIR_00P   ])[k   ]
-			(D.f[DIR_00M])[kb] = mfbbc;//(D.f[ DIR_00M   ])[kb   ] = mfbbc;// -  c2over27 ;	 (D.f[ DIR_00M   ])[kb  ]
-			(D.f[DIR_PP0])[k] = mfaab;//(D.f[ DIR_PP0  ])[kne  ] = mfaab;// -  c1over54 ;	 (D.f[ DIR_PP0  ])[k   ]
-			(D.f[DIR_MM0])[ksw] = mfccb;//(D.f[ DIR_MM0  ])[ksw  ] = mfccb;// -  c1over54 ;	 (D.f[ DIR_MM0  ])[ksw ]
-			(D.f[DIR_PM0])[ks] = mfacb;//(D.f[ DIR_PM0  ])[kse  ] = mfacb;// -  c1over54 ;	 (D.f[ DIR_PM0  ])[ks  ]
-			(D.f[DIR_MP0])[kw] = mfcab;//(D.f[ DIR_MP0  ])[knw  ] = mfcab;// -  c1over54 ;	 (D.f[ DIR_MP0  ])[kw  ]
-			(D.f[DIR_P0P])[k] = mfaba;//(D.f[ DIR_P0P  ])[kte  ] = mfaba;// -  c1over54 ;	 (D.f[ DIR_P0P  ])[k   ]
-			(D.f[DIR_M0M])[kbw] = mfcbc;//(D.f[ DIR_M0M  ])[kbw  ] = mfcbc;// -  c1over54 ;	 (D.f[ DIR_M0M  ])[kbw ]
-			(D.f[DIR_P0M])[kb] = mfabc;//(D.f[ DIR_P0M  ])[kbe  ] = mfabc;// -  c1over54 ;	 (D.f[ DIR_P0M  ])[kb  ]
-			(D.f[DIR_M0P])[kw] = mfcba;//(D.f[ DIR_M0P  ])[ktw  ] = mfcba;// -  c1over54 ;	 (D.f[ DIR_M0P  ])[kw  ]
-			(D.f[DIR_0PP])[k] = mfbaa;//(D.f[ DIR_0PP  ])[ktn  ] = mfbaa;// -  c1over54 ;	 (D.f[ DIR_0PP  ])[k   ]
-			(D.f[DIR_0MM])[kbs] = mfbcc;//(D.f[ DIR_0MM  ])[kbs  ] = mfbcc;// -  c1over54 ;	 (D.f[ DIR_0MM  ])[kbs ]
-			(D.f[DIR_0PM])[kb] = mfbac;//(D.f[ DIR_0PM  ])[kbn  ] = mfbac;// -  c1over54 ;	 (D.f[ DIR_0PM  ])[kb  ]
-			(D.f[DIR_0MP])[ks] = mfbca;//(D.f[ DIR_0MP  ])[kts  ] = mfbca;// -  c1over54 ;	 (D.f[ DIR_0MP  ])[ks  ]
-			(D.f[DIR_000])[k] = mfbbb;//(D.f[ DIR_000])[kzero] = mfbbb;// -  c8over27 ;	 (D.f[ DIR_000])[k   ]
-			(D.f[DIR_PPP])[k] = mfaaa;//(D.f[ DIR_PPP ])[ktne ] = mfaaa;// -  c1over216;	 (D.f[ DIR_PPP ])[k   ]
-			(D.f[DIR_PMP])[ks] = mfaca;//(D.f[ DIR_PMP ])[ktse ] = mfaca;// -  c1over216;	 (D.f[ DIR_PMP ])[ks  ]
-			(D.f[DIR_PPM])[kb] = mfaac;//(D.f[ DIR_PPM ])[kbne ] = mfaac;// -  c1over216;	 (D.f[ DIR_PPM ])[kb  ]
-			(D.f[DIR_PMM])[kbs] = mfacc;//(D.f[ DIR_PMM ])[kbse ] = mfacc;// -  c1over216;	 (D.f[ DIR_PMM ])[kbs ]
-			(D.f[DIR_MPP])[kw] = mfcaa;//(D.f[ DIR_MPP ])[ktnw ] = mfcaa;// -  c1over216;	 (D.f[ DIR_MPP ])[kw  ]
-			(D.f[DIR_MMP])[ksw] = mfcca;//(D.f[ DIR_MMP ])[ktsw ] = mfcca;// -  c1over216;	 (D.f[ DIR_MMP ])[ksw ]
-			(D.f[DIR_MPM])[kbw] = mfcac;//(D.f[ DIR_MPM ])[kbnw ] = mfcac;// -  c1over216;	 (D.f[ DIR_MPM ])[kbw ]
-			(D.f[DIR_MMM])[kbsw] = mfccc;//(D.f[ DIR_MMM ])[kbsw ] = mfccc;// -  c1over216;	 (D.f[ DIR_MMM ])[kbsw]
-										////////////////////////////////////////////////////////////////////////////////////
-		}
-	}
-}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp_Device.cuh b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp_Device.cuh
deleted file mode 100644
index f44842057d554498b0b5d4c733e2425e524a3b75..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp_Device.cuh
+++ /dev/null
@@ -1,18 +0,0 @@
-#ifndef LB_Kernel_CUMULANT_K17_COMP_H
-#define LB_Kernel_CUMULANT_K17_COMP_H
-
-#include <DataTypes.h>
-#include <curand.h>
-
-__global__ void LB_Kernel_CumulantK17Comp(	real omega,
-														unsigned int* bcMatD,
-														unsigned int* neighborX,
-														unsigned int* neighborY,
-														unsigned int* neighborZ,
-														real* DDStart,
-														int size_Mat,
-														int level,
-														real* forces,
-                                                        real* quadricLimiters,
-														bool EvenOrOdd);
-#endif
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty_Device.cu b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17_Device.cu
similarity index 82%
rename from src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty_Device.cu
rename to src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17_Device.cu
index 92aff766b9f8675cc31ce85464cf1848e54bcade..d2b679395a673b5b1d74e3f499a1d53b154d4b89 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty_Device.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17_Device.cu
@@ -1,4 +1,4 @@
-//=======================================================================================
+
 // ____          ____    __    ______     __________   __      __       __        __         
 // \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |        
 //  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |        
@@ -26,17 +26,20 @@
 //  You should have received a copy of the GNU General Public License along
 //  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
 //
-//! \file CumulantK17Almighty_Device.cu
-//! \author Henry Korb, Henrik Asmuth
-//! \date 16/05/2022
-//! \brief CumulantK17CompChim kernel by Martin Schönherr that inlcudes turbulent viscosity and other small mods.
-//!
-//! Additions to CumulantK17CompChim:
-//!     - can incorporate local body force 
-//!     - when applying a local body force, the total round of error of forcing+bodyforce is saved and added in next time step
-//!     - uses turbulent viscosity that is computed in separate kernel (as of now AMD)
-//!     - saves macroscopic values (needed for instance for probes, AMD, and actuator models)
+//! \file CumlantK17_Device.cu
+//! \author Anna Wellmann, Martin Schönherr, Henry Korb, Henrik Asmuth
+//! \date 05/12/2022
+//! \brief Kernel for CumulantK17 including different turbulence models and options for local body forces and writing macroscopic variables
 //!
+//! CumulantK17 kernel using chimera transformations and quartic limiters as present in Geier et al. (2017). Additional options are three different
+//! eddy-viscosity turbulence models (Smagorinsky, AMD, QR) that can be set via the template parameter turbulenceModel (with default 
+//! TurbulenceModel::None). 
+//! The kernel is executed separately for each subset of fluid node indices with a different tag CollisionTemplate. For each subset, only the locally
+//! required options are switched on ( \param writeMacroscopicVariables and/or \param applyBodyForce) in order to minimize memory accesses. The default 
+//! refers to the plain cumlant kernel (CollisionTemplate::Default). 
+//! Nodes are added to subsets (taggedFluidNodes) in Simulation::init using a corresponding tag with different values of CollisionTemplate. These subsets 
+//! are provided by the utilized PostCollisionInteractiors depending on they specifc requirements (e.g. writeMacroscopicVariables for probes).
+
 //=======================================================================================
 /* Device code */
 #include "LBM/LB.h" 
@@ -53,7 +56,7 @@ using namespace vf::lbm::dir;
 
 ////////////////////////////////////////////////////////////////////////////////
 template<TurbulenceModel turbulenceModel, bool writeMacroscopicVariables, bool applyBodyForce>
-__global__ void LB_Kernel_CumulantK17Almighty(
+__global__ void LB_Kernel_CumulantK17(
 	real omega_in,
 	uint* typeOfGridNode,
 	uint* neighborX,
@@ -232,7 +235,7 @@ __global__ void LB_Kernel_CumulantK17Almighty(
         vvy += acc_y;
         vvz += acc_z;
         
-        //    // Reset body force. To be used when not using round-off correction.
+        // Reset body force. To be used when not using round-off correction.
         bodyForceX[k_000] = 0.0f;
         bodyForceY[k_000] = 0.0f;
         bodyForceZ[k_000] = 0.0f;
@@ -692,34 +695,34 @@ __global__ void LB_Kernel_CumulantK17Almighty(
     (dist.f[DIR_MMM])[k_MMM]  = f_PPP;
 }
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::AMD, true, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::AMD, true, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::Smagorinsky, true, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::Smagorinsky, true, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::QR, true, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::QR, true, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::None, true, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::None, true, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::AMD, true, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::AMD, true, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::Smagorinsky, true, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::Smagorinsky, true, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::QR, true, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::QR, true, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::None, true, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::None, true, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::AMD, false, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::AMD, false, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::Smagorinsky, false, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::Smagorinsky, false, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::QR, false, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::QR, false, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::None, false, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::None, false, true > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::AMD, false, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::AMD, false, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::Smagorinsky, false, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::Smagorinsky, false, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::QR, false, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::QR, false, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
 
-template __global__ void LB_Kernel_CumulantK17Almighty < TurbulenceModel::None, false, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
\ No newline at end of file
+template __global__ void LB_Kernel_CumulantK17 < TurbulenceModel::None, false, false > ( 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, const uint *fluidNodeIndices, uint numberOfFluidNodes);
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty_Device.cuh b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17_Device.cuh
similarity index 86%
rename from src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty_Device.cuh
rename to src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17_Device.cuh
index b70b2a75669e1dc12eb3bbd9e1f150fd9c636d63..9d56098f5053cc765b12eab0244a890d18209b1b 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty_Device.cuh
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17_Device.cuh
@@ -1,10 +1,10 @@
-#ifndef LB_Kernel_CUMULANT_K17_ALMIGHTY_H
-#define LB_Kernel_CUMULANT_K17_ALMIGHTY_H
+#ifndef LB_Kernel_CUMULANT_K17_H
+#define LB_Kernel_CUMULANT_K17_H
 
 #include <DataTypes.h>
 #include <curand.h>
 
-template< TurbulenceModel turbulenceModel, bool writeMacroscopicVariables, bool applyBodyForce > __global__ void LB_Kernel_CumulantK17Almighty(
+template< TurbulenceModel turbulenceModel, bool writeMacroscopicVariables, bool applyBodyForce > __global__ void LB_Kernel_CumulantK17(
 	real omega_in,
 	uint* typeOfGridNode,
 	uint* neighborX,
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty.h b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty.h
deleted file mode 100644
index 1819537c095c0106b707ab551b798c0ed62c0e41..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty.h
+++ /dev/null
@@ -1,20 +0,0 @@
-#ifndef CUMULANT_K17_ALMIGHTY_H
-#define CUMULANT_K17_ALMIGHTY_H
-
-#include "Kernel/KernelImp.h"
-#include "Parameter/Parameter.h"
-
-template<TurbulenceModel turbulenceModel> 
-class CumulantK17Almighty : public KernelImp
-{
-public:
-	static std::shared_ptr< CumulantK17Almighty<turbulenceModel> > getNewInstance(std::shared_ptr< Parameter> para, int level);
-	void run() override;
-    void runOnIndices(const unsigned int *indices, unsigned int size_indices, CollisionTemplate collisionTemplate, CudaStreamIndex streamIndex) override;
-
-private:
-    CumulantK17Almighty();
-    CumulantK17Almighty(std::shared_ptr<Parameter> para, int level);
-};
-
-#endif 
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.cpp b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.cpp
index 553b2b508527f3f94beeeff71334124410a16864..5a2d8c9a426e5cb23ca75f91aaf6fbff75cba72b 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.cpp
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.cpp
@@ -8,10 +8,9 @@
 #include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/BGKPlus/BGKPlusCompSP27.h"
 #include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/Cascade/CascadeCompSP27.h"
 #include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/Cumulant/CumulantCompSP27.h"
-#include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Comp.h"
 #include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Unified/CumulantK17Unified.h"
 #include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17chim/CumulantK17CompChim.h"
-#include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Almighty/CumulantK17Almighty.h"
+#include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17.h"
 #include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17Bulk/CumulantK17BulkComp.h"
 #include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantAll4/CumulantAll4CompSP27.h"
 #include "Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK18/CumulantK18Comp.h"
@@ -113,9 +112,6 @@ std::shared_ptr<Kernel> KernelFactoryImp::makeKernel(std::shared_ptr<Parameter>
     } else if (kernel == "CumulantCompSP27") {
         newKernel     = CumulantCompSP27::getNewInstance(para, level);
         checkStrategy = FluidFlowCompStrategy::getInstance();
-    } else if (kernel == "CumulantK17Comp") {
-        newKernel     = CumulantK17Comp::getNewInstance(para, level);
-        checkStrategy = FluidFlowCompStrategy::getInstance();
     } else if (kernel == "CumulantK15Unified") {
         newKernel     = std::make_shared<vf::gpu::CumulantK15Unified>(para, level);
         checkStrategy = FluidFlowCompStrategy::getInstance();
@@ -128,20 +124,20 @@ std::shared_ptr<Kernel> KernelFactoryImp::makeKernel(std::shared_ptr<Parameter>
     } else if (kernel == "CumulantK17CompChim") {
         newKernel     = CumulantK17CompChim::getNewInstance(para, level);
         checkStrategy = FluidFlowCompStrategy::getInstance();
-    } else if (kernel == "CumulantK17Almighty"){               
+    } else if (kernel == "CumulantK17"){               
         switch(para->getTurbulenceModel())                                          
         {   
             case TurbulenceModel::AMD:
-                newKernel = CumulantK17Almighty<TurbulenceModel::AMD>::getNewInstance(para, level);   
+                newKernel = CumulantK17<TurbulenceModel::AMD>::getNewInstance(para, level);   
                 break;
             case TurbulenceModel::Smagorinsky:
-                newKernel = CumulantK17Almighty<TurbulenceModel::Smagorinsky>::getNewInstance(para, level);  
+                newKernel = CumulantK17<TurbulenceModel::Smagorinsky>::getNewInstance(para, level);  
                 break;
             case TurbulenceModel::QR:
-                newKernel = CumulantK17Almighty<TurbulenceModel::QR>::getNewInstance(para, level);  
+                newKernel = CumulantK17<TurbulenceModel::QR>::getNewInstance(para, level);  
                 break;
             case TurbulenceModel::None:
-                newKernel = CumulantK17Almighty<TurbulenceModel::None>::getNewInstance(para, level); 
+                newKernel = CumulantK17<TurbulenceModel::None>::getNewInstance(para, level); 
                 break;
             default:
                 throw std::runtime_error("Unknown turbulence model!");
diff --git a/src/gpu/VirtualFluids_GPU/LBM/LB.h b/src/gpu/VirtualFluids_GPU/LBM/LB.h
index bd9ce06a01bddf5504ed15267801b409d0a08461..904471123a895d65f33c8d91e6c5e5ed0296a9f6 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/LB.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/LB.h
@@ -64,7 +64,7 @@ enum class TurbulenceModel {
    None
 };
 
-//! \brief An enumeration for selecting a template of the collision kernel (CumulantK17Almighty)
+//! \brief An enumeration for selecting a template of the collision kernel (CumulantK17)
 enum class CollisionTemplate {
    //! - Default: plain collision without additional read/write
    Default,
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index 8ff0fdf85fe4aab2cc106fe05a32d81cab28c8d7..b120dd20531e2a4f672f16f59b4e2dfe6decaaf1 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -125,6 +125,7 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     //////////////////////////////////////////////////////////////////////////
     allocNeighborsOffsetsScalesAndBoundaries(gridProvider);
 
+    //! Get tagged fluid nodes with corresponding value for CollisionTemplate from interactors
     for (SPtr<PreCollisionInteractor> actuator : para->getActuators()) {
         actuator->init(para.get(), &gridProvider, cudaMemoryManager.get());
         actuator->getTaggedFluidNodes( para.get(), &gridProvider );
@@ -144,6 +145,7 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
         para->getStreamManager()->createCudaEvents();
     }
     //////////////////////////////////////////////////////////////////////////
+    
     if (para->getKernelNeedsFluidNodeIndicesToRun())
     {
         gridProvider.sortFluidNodeTags();
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
index 6505dd32399726bf01e54446dc22b03b37578d6e..bcc3c423f36407a78d1120ee64c03ab17d8b9a6d 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -1613,7 +1613,7 @@ void Parameter::setOutflowBoundaryNormalZ(std::string outflowNormalZ)
 void Parameter::setMainKernel(std::string kernel)
 {
     this->mainKernel = kernel;
-    if ( kernel.find("Almighty") != std::string::npos )
+    if ( kernel.find("CumulantK17") != std::string::npos )
         this->kernelNeedsFluidNodeIndicesToRun = true;
 }
 void Parameter::setMultiKernelOn(bool isOn)
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/PrecursorWriter.h b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/PrecursorWriter.h
index 03df46b8a8437e5e9e82f1e28c2d833ee424430a..8dc8e36fe57713513c58ed76637ab3fdc34a251f 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/PrecursorWriter.h
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/PrecursorWriter.h
@@ -1,3 +1,41 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __         
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |        
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |        
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |        
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____    
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|   
+//      \    \  |    |   ________________________________________________________________    
+//       \    \ |    |  |  ______________________________________________________________|   
+//        \    \|    |  |  |         __          __     __     __     ______      _______    
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)   
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______    
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/   
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file PrecursorWriter.h
+//! \author Henry Korb, Henrik Asmuth
+//! \date 05/12/2022
+//! \brief Probe writing planes of data to be used as inflow data in successor simulation using PrecursorBC
+//!
+//! The probe writes out yz-planes at a specifc x position ( \param xPos ) of either velocity or distributions 
+//! that can be read by PrecursorBC as inflow data.
+//=======================================================================================
+
+
 #ifndef PRECURSORPROBE_H_
 #define PRECURSORPROBE_H_