diff --git a/src/VirtualFluids_GPU/Calculation/ForceCalculations.cpp b/src/VirtualFluids_GPU/Calculation/ForceCalculations.cpp
index 323bebcef462e2f971e83504622a209a0444664f..0f39435b15ca8ebab4e93db14f1bf802532da317 100644
--- a/src/VirtualFluids_GPU/Calculation/ForceCalculations.cpp
+++ b/src/VirtualFluids_GPU/Calculation/ForceCalculations.cpp
@@ -17,8 +17,8 @@ ForceCalculations::ForceCalculations(Parameter* para)
 
 	Kpcrit = 3.0 / Ta;// 0.3;
 	Tcrit = 3.0 * Ta; // 30.0;
-	Tn = 0.5 * Tcrit; //1.0 * Tcrit; //0.5 * Tcrit;
-	Tv = 0.12 * Tcrit; //0.24 * Tcrit; //0.12 * Tcrit;
+	Tn = 1.0 * Tcrit; //0.5 * Tcrit;
+	Tv = 0.24 * Tcrit; //0.12 * Tcrit;
 
 	Kp = 0.6 * Kpcrit;
 	Ki = Kp / Tn;
diff --git a/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp b/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..83e7eb0ae36a758bca75e28a675a0de11881aaa1
--- /dev/null
+++ b/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp
@@ -0,0 +1,103 @@
+#include "Calculation/PorousMedia.h"
+
+//////////////////////////////////////////////////////////////////////////
+//#include "GPU/GPU_Interface.h"
+//#include <cuda_runtime.h>
+//#include <helper_cuda.h>
+//#include <stdio.h>
+//#include <fstream>
+//#include <sstream>
+//using namespace std;
+//////////////////////////////////////////////////////////////////////////
+
+PorousMedia::PorousMedia()
+{
+	porosity = 0.0;
+	geoID = 0;
+	setLBMvaluesToZero();
+	setSIvaluesToZero();
+	setCoordinatesToZero();
+}
+
+PorousMedia::PorousMedia(double porosity, unsigned int geoID, double darcySI, double forchheimerSI, double dxLBM, double dtLBM)
+{
+	this->porosity = porosity;
+	this->geoID = geoID;
+	this->darcySI = darcySI;
+	this->forchheimerSI = forchheimerSI;
+	this->dxLBM = dxLBM;
+	this->dtLBM = dtLBM;
+	this->forchheimerLBM = this->forchheimerSI * this->dtLBM;
+	this->darcyLBM = this->darcySI * this->dxLBM * this->dxLBM;
+	setCoordinatesToZero();
+}
+
+PorousMedia::~PorousMedia()
+{
+}
+
+
+//setter
+void PorousMedia::setPorosity(double porosity) {			this->porosity = porosity; }
+void PorousMedia::setDarcySI(double darcySI) {				this->darcySI = darcySI; }
+void PorousMedia::setForchheimerSI(double forchheimerSI) {	this->forchheimerSI = forchheimerSI; }
+
+void PorousMedia::setStartCoordinates(double startCoordX, double startCoordY, double startCoordZ) 
+{
+	this->startCoordX = startCoordX;
+	this->startCoordY = startCoordY;
+	this->startCoordZ = startCoordZ;
+}
+
+void PorousMedia::setEndCoordinates(double endCoordX, double endCoordY, double endCoordZ)
+{
+	this->endCoordX = endCoordX;
+	this->endCoordY = endCoordY;
+	this->endCoordZ = endCoordZ;
+}
+
+void PorousMedia::setHostNodeIDsPM(unsigned int* hostNodeIDsPM)
+{
+	this->hostNodeIDsPM = hostNodeIDsPM;
+}
+
+void PorousMedia::setDeviceNodeIDsPM(unsigned int* deviceNodeIDsPM)
+{
+	this->deviceNodeIDsPM = deviceNodeIDsPM;
+}
+
+//void PorousMedia::definePMarea(Parameter* para, unsigned int level)
+//{
+//	unsigned int counter = 0;
+//	for (unsigned int i = 0; i < para->getParH(level)->size_Mat_SP; i++)
+//	{
+//		if (((para->getParH(level)->coordX_SP[i] >= this->startCoordX) && (para->getParH(level)->coordX_SP[i] <= this->endCoordX)) &&
+//			((para->getParH(level)->coordY_SP[i] >= this->startCoordY) && (para->getParH(level)->coordY_SP[i] <= this->endCoordY)) && 
+//			((para->getParH(level)->coordZ_SP[i] >= this->startCoordZ) && (para->getParH(level)->coordZ_SP[i] <= this->endCoordZ)) )
+//		{
+//			if (para->getParH(level)->geoSP[i] >= GEO_FLUID)
+//			{
+//				para->getParH(level)->geoSP[i] = this->geoID;
+//				nodeIDsPorousMedia.push_back(i);
+//				counter++;
+//			}
+//		}
+//	}
+//
+//	this->sizePM = counter;
+//
+//	for (unsigned int j = 0; j <= this->sizePM; j++)
+//	{
+//	}
+//}
+
+//getter
+double PorousMedia::getPorosity(){					return this->porosity; }
+double PorousMedia::getDarcySI() {					return this->darcySI; }
+double PorousMedia::getForchheimerSI(){				return this->forchheimerSI; }
+double PorousMedia::getDarcyLBM() {					return this->darcyLBM; }
+double PorousMedia::getForchheimerLBM() {			return this->forchheimerLBM; }
+unsigned int PorousMedia::getGeoID() {				return this->geoID; }
+unsigned int PorousMedia::getSizePM() {				return this->sizePM; }
+unsigned int* PorousMedia::getHostNodeIDsPM() {		return this->hostNodeIDsPM; }
+unsigned int* PorousMedia::getDeviceNodeIDsPM() {	return this->deviceNodeIDsPM; }
diff --git a/src/VirtualFluids_GPU/Calculation/PorousMedia.h b/src/VirtualFluids_GPU/Calculation/PorousMedia.h
new file mode 100644
index 0000000000000000000000000000000000000000..da0c9741b7501affb6284e8609157afc811d808b
--- /dev/null
+++ b/src/VirtualFluids_GPU/Calculation/PorousMedia.h
@@ -0,0 +1,89 @@
+#ifndef POROUS_MEDIA_H
+#define POROUS_MEDIA_H
+
+//#include "Parameter/Parameter.h"
+//#include "Utilities/StringUtil.hpp"
+//#include "basics/utilities/UbSystem.h"
+
+#include <iostream>
+#include <vector>
+#include <stdio.h>
+using namespace std;
+
+class PorousMedia
+{
+public:
+	PorousMedia();
+	PorousMedia(double porosity, unsigned int geoID, double darcySI, double forchheimerSI, double dxLBM, double dtLBM);
+	~PorousMedia();
+
+	//setter
+	void setPorosity(double porosity);
+	void setDarcySI(double darcySI);
+	void setForchheimerSI(double forchheimerSI);
+	void setStartCoordinates(double startCoordX, double startCoordY, double startCoordZ);
+	void setEndCoordinates(double endCoordX, double endCoordY, double endCoordZ);
+	//void definePMarea(Parameter* para, unsigned int level);
+	void setHostNodeIDsPM(unsigned int* hostNodeIDsPM);
+	void setDeviceNodeIDsPM(unsigned int* deviceNodeIDsPM);
+
+	//getter
+	double getPorosity();
+	double getDarcySI();
+	double getForchheimerSI();
+	double getDarcyLBM();
+	double getForchheimerLBM();
+	unsigned int getGeoID();
+	unsigned int getSizePM();
+	unsigned int* getHostNodeIDsPM();
+	unsigned int* getDeviceNodeIDsPM();
+
+private:
+	double porosity;
+	double darcySI; //[1/s]
+	double darcyLBM; 
+	double forchheimerSI; //[1/m]
+	double forchheimerLBM;
+	double dxLBM;
+	double dtLBM;
+	double startCoordX;
+	double startCoordY;
+	double startCoordZ;
+	double endCoordX;
+	double endCoordY;
+	double endCoordZ;
+	unsigned int geoID;
+	std::vector< unsigned int > nodeIDsPorousMedia;
+	unsigned int sizePM;
+	unsigned int *hostNodeIDsPM;
+	unsigned int *deviceNodeIDsPM;
+
+	void setCoordinatesToZero() 
+	{
+		startCoordX = 0;
+		startCoordY = 0;
+		startCoordZ = 0;
+		endCoordX = 0;
+		endCoordY = 0;
+		endCoordZ = 0;
+	};
+
+	void setSIvaluesToZero()
+	{
+		darcySI = 0.0;
+		forchheimerSI = 0.0;
+	};
+
+	void setLBMvaluesToZero()
+	{
+		darcyLBM = 0.0;
+		forchheimerLBM = 0.0;
+		dxLBM = 0.0;
+		dtLBM = 0.0;
+	};
+
+
+};
+
+
+#endif /* POROUS_MEDIA_H */
diff --git a/src/VirtualFluids_GPU/GPU/GPU_Interface.h b/src/VirtualFluids_GPU/GPU/GPU_Interface.h
index b882422e1875ed7a20049979798bfa8dbb2b0535..e3cd75a7587b2d70ddbfbca77a653947d339f9bf 100644
--- a/src/VirtualFluids_GPU/GPU/GPU_Interface.h
+++ b/src/VirtualFluids_GPU/GPU/GPU_Interface.h
@@ -311,6 +311,22 @@ extern "C" void KernelWaleCumAA2016CompSP27( unsigned int numberOfThreads,
 											 doubflo* forces,
 											 bool EvenOrOdd);
 
+extern "C" void KernelPMCumOneCompSP27(unsigned int numberOfThreads, 
+									   doubflo omega,
+									   unsigned int* bcMatD,
+									   unsigned int* neighborX,
+									   unsigned int* neighborY,
+									   unsigned int* neighborZ,
+									   doubflo* DD,
+									   int size_Mat,
+									   int level,
+									   doubflo* forces,
+									   doubflo porosity,
+									   doubflo darcy,
+									   doubflo forchheimer,
+									   unsigned int porousMedia,
+									   bool EvenOrOdd);
+
 extern "C" void KernelThS7(unsigned int numberOfThreads, 
                            doubflo diffusivity,
                            unsigned int* bcMatD,
diff --git a/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh b/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
index bc962d489660f28e7112daac152efbc2e10bc61f..4eadea9766a24daa934fc5949d06fc6482cab778 100644
--- a/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
+++ b/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
@@ -278,6 +278,21 @@ extern "C" __global__ void LB_Kernel_Wale_Cum_AA2016_Comp_SP_27( doubflo omega,
 																 doubflo* forces,
 																 bool EvenOrOdd);
 
+extern "C" __global__ void LB_Kernel_PM_Cum_One_Comp_SP_27( doubflo omega,
+															unsigned int* bcMatD,
+															unsigned int* neighborX,
+															unsigned int* neighborY,
+															unsigned int* neighborZ,
+															doubflo* DDStart,
+															int size_Mat,
+															int level,
+															doubflo* forces,
+															doubflo porosity,
+															doubflo darcy,
+															doubflo forchheimer,
+															unsigned int porousMedia,
+															bool EvenOrOdd);
+
 extern "C" __global__ void LB_Kernel_ThS7(doubflo diffusivity,
                                           unsigned int* bcMatD,
                                           unsigned int* neighborX,
diff --git a/src/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/VirtualFluids_GPU/GPU/LBMKernel.cu
index f377dff6f93ed88bd78cc76be2f2097fc4bd84b1..6c4d91e169b37d2705ea3e321c29ca1101bccc20 100644
--- a/src/VirtualFluids_GPU/GPU/LBMKernel.cu
+++ b/src/VirtualFluids_GPU/GPU/LBMKernel.cu
@@ -957,6 +957,54 @@ extern "C" void KernelWaleCumOneCompSP27(unsigned int numberOfThreads,
 		getLastCudaError("LB_Kernel_Wale_Cum_One_Comp_SP_27 execution failed"); 
 }
 //////////////////////////////////////////////////////////////////////////
+extern "C" void KernelPMCumOneCompSP27(unsigned int numberOfThreads, 
+									   doubflo omega,
+									   unsigned int* bcMatD,
+									   unsigned int* neighborX,
+									   unsigned int* neighborY,
+									   unsigned int* neighborZ,
+									   doubflo* DD,
+									   int size_Mat,
+									   int level,
+									   doubflo* forces,
+									   doubflo porosity,
+									   doubflo darcy,
+									   doubflo forchheimer,
+									   unsigned int porousMedia,
+									   bool EvenOrOdd)
+{
+	int Grid = (size_Mat / numberOfThreads) + 1;
+	int Grid1, Grid2;
+	if (Grid > 512)
+	{
+		Grid1 = 512;
+		Grid2 = (Grid / Grid1) + 1;
+	}
+	else
+	{
+		Grid1 = 1;
+		Grid2 = Grid;
+	}
+	dim3 grid(Grid1, Grid2, 1);
+	dim3 threads(numberOfThreads, 1, 1);
+
+	LB_Kernel_PM_Cum_One_Comp_SP_27 <<< grid, threads >>>(omega,
+														  bcMatD,
+														  neighborX,
+														  neighborY,
+														  neighborZ,
+														  DD,
+														  size_Mat,
+														  level,
+														  forces,
+														  porosity,
+														  darcy,
+														  forchheimer,
+														  porousMedia,
+														  EvenOrOdd); 
+	getLastCudaError("LB_Kernel_PM_Cum_One_Comp_SP_27 execution failed"); 
+}
+//////////////////////////////////////////////////////////////////////////
 extern "C" void KernelWaleCumAA2016CompSP27( unsigned int numberOfThreads,
 											 doubflo s9,
 											 unsigned int* bcMatD,
diff --git a/src/VirtualFluids_GPU/GPU/PorousMediaCumulant27.cu b/src/VirtualFluids_GPU/GPU/PorousMediaCumulant27.cu
new file mode 100644
index 0000000000000000000000000000000000000000..fdbc535bbad930a4d2f403e8e725e08eeec42c36
--- /dev/null
+++ b/src/VirtualFluids_GPU/GPU/PorousMediaCumulant27.cu
@@ -0,0 +1,942 @@
+/* Device code */
+#include "LBM/D3Q27.h"
+#include "math.h"
+#include "GPU/constant.h"
+
+////////////////////////////////////////////////////////////////////////////////
+extern "C" __global__ void LB_Kernel_PM_Cum_One_Comp_SP_27( doubflo omega,
+															unsigned int* bcMatD,
+															unsigned int* neighborX,
+															unsigned int* neighborY,
+															unsigned int* neighborZ,
+															doubflo* DDStart,
+															int size_Mat,
+															int level,
+															doubflo* forces,
+															doubflo porosity,
+															doubflo darcy,
+															doubflo forchheimer,
+															unsigned int porousMedia,
+															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 == porousMedia)
+		{
+			Distributions27 D;
+			if (EvenOrOdd==true)
+			{
+				D.f[dirE   ] = &DDStart[dirE   *size_Mat];
+				D.f[dirW   ] = &DDStart[dirW   *size_Mat];
+				D.f[dirN   ] = &DDStart[dirN   *size_Mat];
+				D.f[dirS   ] = &DDStart[dirS   *size_Mat];
+				D.f[dirT   ] = &DDStart[dirT   *size_Mat];
+				D.f[dirB   ] = &DDStart[dirB   *size_Mat];
+				D.f[dirNE  ] = &DDStart[dirNE  *size_Mat];
+				D.f[dirSW  ] = &DDStart[dirSW  *size_Mat];
+				D.f[dirSE  ] = &DDStart[dirSE  *size_Mat];
+				D.f[dirNW  ] = &DDStart[dirNW  *size_Mat];
+				D.f[dirTE  ] = &DDStart[dirTE  *size_Mat];
+				D.f[dirBW  ] = &DDStart[dirBW  *size_Mat];
+				D.f[dirBE  ] = &DDStart[dirBE  *size_Mat];
+				D.f[dirTW  ] = &DDStart[dirTW  *size_Mat];
+				D.f[dirTN  ] = &DDStart[dirTN  *size_Mat];
+				D.f[dirBS  ] = &DDStart[dirBS  *size_Mat];
+				D.f[dirBN  ] = &DDStart[dirBN  *size_Mat];
+				D.f[dirTS  ] = &DDStart[dirTS  *size_Mat];
+				D.f[dirZERO] = &DDStart[dirZERO*size_Mat];
+				D.f[dirTNE ] = &DDStart[dirTNE *size_Mat];
+				D.f[dirTSW ] = &DDStart[dirTSW *size_Mat];
+				D.f[dirTSE ] = &DDStart[dirTSE *size_Mat];
+				D.f[dirTNW ] = &DDStart[dirTNW *size_Mat];
+				D.f[dirBNE ] = &DDStart[dirBNE *size_Mat];
+				D.f[dirBSW ] = &DDStart[dirBSW *size_Mat];
+				D.f[dirBSE ] = &DDStart[dirBSE *size_Mat];
+				D.f[dirBNW ] = &DDStart[dirBNW *size_Mat];
+			}
+			else
+			{
+				D.f[dirW   ] = &DDStart[dirE   *size_Mat];
+				D.f[dirE   ] = &DDStart[dirW   *size_Mat];
+				D.f[dirS   ] = &DDStart[dirN   *size_Mat];
+				D.f[dirN   ] = &DDStart[dirS   *size_Mat];
+				D.f[dirB   ] = &DDStart[dirT   *size_Mat];
+				D.f[dirT   ] = &DDStart[dirB   *size_Mat];
+				D.f[dirSW  ] = &DDStart[dirNE  *size_Mat];
+				D.f[dirNE  ] = &DDStart[dirSW  *size_Mat];
+				D.f[dirNW  ] = &DDStart[dirSE  *size_Mat];
+				D.f[dirSE  ] = &DDStart[dirNW  *size_Mat];
+				D.f[dirBW  ] = &DDStart[dirTE  *size_Mat];
+				D.f[dirTE  ] = &DDStart[dirBW  *size_Mat];
+				D.f[dirTW  ] = &DDStart[dirBE  *size_Mat];
+				D.f[dirBE  ] = &DDStart[dirTW  *size_Mat];
+				D.f[dirBS  ] = &DDStart[dirTN  *size_Mat];
+				D.f[dirTN  ] = &DDStart[dirBS  *size_Mat];
+				D.f[dirTS  ] = &DDStart[dirBN  *size_Mat];
+				D.f[dirBN  ] = &DDStart[dirTS  *size_Mat];
+				D.f[dirZERO] = &DDStart[dirZERO*size_Mat];
+				D.f[dirBSW ] = &DDStart[dirTNE *size_Mat];
+				D.f[dirBNE ] = &DDStart[dirTSW *size_Mat];
+				D.f[dirBNW ] = &DDStart[dirTSE *size_Mat];
+				D.f[dirBSE ] = &DDStart[dirTNW *size_Mat];
+				D.f[dirTSW ] = &DDStart[dirBNE *size_Mat];
+				D.f[dirTNE ] = &DDStart[dirBSW *size_Mat];
+				D.f[dirTNW ] = &DDStart[dirBSE *size_Mat];
+				D.f[dirTSE ] = &DDStart[dirBNW *size_Mat];
+			}
+
+			////////////////////////////////////////////////////////////////////////////////
+			//index
+			unsigned int kw   = neighborX[k];
+			unsigned int ks   = neighborY[k];
+			unsigned int kb   = neighborZ[k];
+			unsigned int ksw  = neighborY[kw];
+			unsigned int kbw  = neighborZ[kw];
+			unsigned int kbs  = neighborZ[ks];
+			unsigned int kbsw = neighborZ[ksw];
+			//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			doubflo mfcbb = (D.f[dirE   ])[k  ];
+			doubflo mfabb = (D.f[dirW   ])[kw ];
+			doubflo mfbcb = (D.f[dirN   ])[k  ];
+			doubflo mfbab = (D.f[dirS   ])[ks ];
+			doubflo mfbbc = (D.f[dirT   ])[k  ];
+			doubflo mfbba = (D.f[dirB   ])[kb ];
+			doubflo mfccb = (D.f[dirNE  ])[k  ];
+			doubflo mfaab = (D.f[dirSW  ])[ksw];
+			doubflo mfcab = (D.f[dirSE  ])[ks ];
+			doubflo mfacb = (D.f[dirNW  ])[kw ];
+			doubflo mfcbc = (D.f[dirTE  ])[k  ];
+			doubflo mfaba = (D.f[dirBW  ])[kbw];
+			doubflo mfcba = (D.f[dirBE  ])[kb ];
+			doubflo mfabc = (D.f[dirTW  ])[kw ];
+			doubflo mfbcc = (D.f[dirTN  ])[k  ];
+			doubflo mfbaa = (D.f[dirBS  ])[kbs];
+			doubflo mfbca = (D.f[dirBN  ])[kb ];
+			doubflo mfbac = (D.f[dirTS  ])[ks ];
+			doubflo mfbbb = (D.f[dirZERO])[k  ];
+			doubflo mfccc = (D.f[dirTNE ])[k  ];
+			doubflo mfaac = (D.f[dirTSW ])[ksw];
+			doubflo mfcac = (D.f[dirTSE ])[ks ];
+			doubflo mfacc = (D.f[dirTNW ])[kw ];
+			doubflo mfcca = (D.f[dirBNE ])[kb ];
+			doubflo mfaaa = (D.f[dirBSW])[kbsw];
+			doubflo mfcaa = (D.f[dirBSE ])[kbs];
+			doubflo mfaca = (D.f[dirBNW ])[kbw];
+			////////////////////////////////////////////////////////////////////////////////////
+			doubflo 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;
+
+			doubflo rho = one+drho;
+			////////////////////////////////////////////////////////////////////////////////////
+			doubflo vvx    =((((mfccc-mfaaa) + (mfcac-mfaca)) + ((mfcaa-mfacc) + (mfcca-mfaac))) + 
+						     (((mfcba-mfabc) + (mfcbc-mfaba)) + ((mfcab-mfacb) + (mfccb-mfaab))) +
+						       (mfcbb-mfabb)) / rho;
+			doubflo vvy    =((((mfccc-mfaaa) + (mfaca-mfcac)) + ((mfacc-mfcaa) + (mfcca-mfaac))) + 
+				             (((mfbca-mfbac) + (mfbcc-mfbaa)) + ((mfacb-mfcab) + (mfccb-mfaab))) +
+				               (mfbcb-mfbab)) / rho;
+			doubflo vvz    =((((mfccc-mfaaa) + (mfcac-mfaca)) + ((mfacc-mfcaa) + (mfaac-mfcca))) + 
+				             (((mfbac-mfbca) + (mfbcc-mfbaa)) + ((mfabc-mfcba) + (mfcbc-mfaba))) +
+				               (mfbbc-mfbba)) / rho;
+			doubflo vx2    = vvx*vvx;
+			doubflo vy2    = vvy*vvy;
+			doubflo vz2    = vvz*vvz;
+			////////////////////////////////////////////////////////////////////////////////////
+			//porous media
+			vvx = -(two * vvx) / (-two + darcy + forchheimer * sqrtf(vx2 + vy2 + vz2));
+			vvy = -(two * vvy) / (-two + darcy + forchheimer * sqrtf(vx2 + vy2 + vz2));
+			vvz = -(two * vvz) / (-two + darcy + forchheimer * sqrtf(vx2 + vy2 + vz2));
+			//vvx = (two * vvx) / (two + 134.4 + 0.0068287 * sqrtf(vx2 + vy2 + vz2));
+			//vvy = (two * vvy) / (two + 134.4 + 0.0068287 * sqrtf(vx2 + vy2 + vz2));
+			//vvz = (two * vvz) / (two + 134.4 + 0.0068287 * sqrtf(vx2 + vy2 + vz2));
+			////////////////////////////////////////////////////////////////////////////////////
+			//the force be with you
+			doubflo fx = forces[0]/(pow(two,level)); //zero;//0.0032653/(pow(two,level)); //0.000000005;//(two/1600000.0) / 120.0; //
+			doubflo fy = forces[1]/(pow(two,level)); //zero;
+			doubflo fz = forces[2]/(pow(two,level)); //zero;
+			vvx += fx*c1o2;
+			vvy += fy*c1o2;
+			vvz += fz*c1o2;
+			////////////////////////////////////////////////////////////////////////////////////
+			//doubflo omega = omega_in;
+			////////////////////////////////////////////////////////////////////////////////////
+			doubflo oMdrho = one; // comp special
+			////////////////////////////////////////////////////////////////////////////////////
+			doubflo m0, m1, m2;	
+			//////////////////////////////////////////////////////////////////////////////////////
+			doubflo wadjust;
+			doubflo qudricLimitP = 0.01f;// * 0.0001f;
+			doubflo qudricLimitM = 0.01f;// * 0.0001f;
+			doubflo qudricLimitD = 0.01f;// * 0.001f;
+			////////////////////////////////////////////////////////////////////////////////////
+			//Hin
+			////////////////////////////////////////////////////////////////////////////////////
+			// mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36  Konditionieren
+			////////////////////////////////////////////////////////////////////////////////////
+			// Z - Dir
+			m2    = mfaaa	+ mfaac;
+			m1    = mfaac	- mfaaa;
+			m0    = m2		+ mfaab;
+			mfaaa = m0;
+			m0   += c1o36 * oMdrho;	
+			mfaab = m1 -		m0 * vvz;
+			mfaac = m2 - two*	m1 * vvz + vz2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfaba  + mfabc;
+			m1    = mfabc  - mfaba;
+			m0    = m2		+ mfabb;
+			mfaba = m0;
+			m0   += c1o9 * oMdrho;
+			mfabb = m1 -		m0 * vvz;
+			mfabc = m2 - two*	m1 * vvz + vz2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfaca  + mfacc;
+			m1    = mfacc  - mfaca;
+			m0    = m2		+ mfacb;
+			mfaca = m0;
+			m0   += c1o36 * oMdrho;
+			mfacb = m1 -		m0 * vvz;
+			mfacc = m2 - two*	m1 * vvz + vz2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfbaa	+ mfbac;
+			m1    = mfbac	- mfbaa;
+			m0    = m2		+ mfbab;
+			mfbaa = m0;
+			m0   += c1o9 * oMdrho;
+			mfbab = m1 -		m0 * vvz;
+			mfbac = m2 - two*	m1 * vvz + vz2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfbba  + mfbbc;
+			m1    = mfbbc  - mfbba;
+			m0    = m2		+ mfbbb;
+			mfbba = m0;
+			m0   += c4o9 * oMdrho;
+			mfbbb = m1 -		m0 * vvz;
+			mfbbc = m2 - two*	m1 * vvz + vz2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfbca  + mfbcc;
+			m1    = mfbcc  - mfbca;
+			m0    = m2		+ mfbcb;
+			mfbca = m0;
+			m0   += c1o9 * oMdrho;
+			mfbcb = m1 -		m0 * vvz;
+			mfbcc = m2 - two*	m1 * vvz + vz2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfcaa	+ mfcac;
+			m1    = mfcac	- mfcaa;
+			m0    = m2		+ mfcab;
+			mfcaa = m0;
+			m0   += c1o36 * oMdrho;
+			mfcab = m1 -		m0 * vvz;
+			mfcac = m2 - two*	m1 * vvz + vz2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfcba  + mfcbc;
+			m1    = mfcbc  - mfcba;
+			m0    = m2		+ mfcbb;
+			mfcba = m0;
+			m0   += c1o9 * oMdrho;
+			mfcbb = m1 -		m0 * vvz;
+			mfcbc = m2 - two*	m1 * vvz + vz2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfcca  + mfccc;
+			m1    = mfccc  - mfcca;
+			m0    = m2		+ mfccb;
+			mfcca = m0;
+			m0   += c1o36 * oMdrho;
+			mfccb = m1 -		m0 * vvz;
+			mfccc = m2 - two*	m1 * vvz + vz2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			// mit  1/6, 0, 1/18, 2/3, 0, 2/9, 1/6, 0, 1/18 Konditionieren
+			////////////////////////////////////////////////////////////////////////////////////
+			// Y - Dir
+			m2    = mfaaa	+ mfaca;
+			m1    = mfaca	- mfaaa;
+			m0    = m2		+ mfaba;
+			mfaaa = m0;
+			m0   += c1o6 * oMdrho;
+			mfaba = m1 -		m0 * vvy;
+			mfaca = m2 - two*	m1 * vvy + vy2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfaab  + mfacb;
+			m1    = mfacb  - mfaab;
+			m0    = m2		+ mfabb;
+			mfaab = m0;
+			mfabb = m1 -		m0 * vvy;
+			mfacb = m2 - two*	m1 * vvy + vy2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfaac  + mfacc;
+			m1    = mfacc  - mfaac;
+			m0    = m2		+ mfabc;
+			mfaac = m0;
+			m0   += c1o18 * oMdrho;
+			mfabc = m1 -		m0 * vvy;
+			mfacc = m2 - two*	m1 * vvy + vy2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfbaa	+ mfbca;
+			m1    = mfbca	- mfbaa;
+			m0    = m2		+ mfbba;
+			mfbaa = m0;
+			m0   += c2o3 * oMdrho;
+			mfbba = m1 -		m0 * vvy;
+			mfbca = m2 - two*	m1 * vvy + vy2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfbab  + mfbcb;
+			m1    = mfbcb  - mfbab;
+			m0    = m2		+ mfbbb;
+			mfbab = m0;
+			mfbbb = m1 -		m0 * vvy;
+			mfbcb = m2 - two*	m1 * vvy + vy2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfbac  + mfbcc;
+			m1    = mfbcc  - mfbac;
+			m0    = m2		+ mfbbc;
+			mfbac = m0;
+			m0   += c2o9 * oMdrho;
+			mfbbc = m1 -		m0 * vvy;
+			mfbcc = m2 - two*	m1 * vvy + vy2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfcaa	+ mfcca;
+			m1    = mfcca	- mfcaa;
+			m0    = m2		+ mfcba;
+			mfcaa = m0;
+			m0   += c1o6 * oMdrho;
+			mfcba = m1 -		m0 * vvy;
+			mfcca = m2 - two*	m1 * vvy + vy2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfcab  + mfccb;
+			m1    = mfccb  - mfcab;
+			m0    = m2		+ mfcbb;
+			mfcab = m0;
+			mfcbb = m1 -		m0 * vvy;
+			mfccb = m2 - two*	m1 * vvy + vy2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfcac  + mfccc;
+			m1    = mfccc  - mfcac;
+			m0    = m2		+ mfcbc;
+			mfcac = m0;
+			m0   += c1o18 * oMdrho;
+			mfcbc = m1 -		m0 * vvy;
+			mfccc = m2 - two*	m1 * vvy + vy2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			// mit     1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9		Konditionieren
+			////////////////////////////////////////////////////////////////////////////////////
+			// X - Dir
+			m2    = mfaaa	+ mfcaa;
+			m1    = mfcaa	- mfaaa;
+			m0    = m2		+ mfbaa;
+			mfaaa = m0;
+			m0   += one* oMdrho;
+			mfbaa = m1 -		m0 * vvx;
+			mfcaa = m2 - two*	m1 * vvx + vx2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfaba  + mfcba;
+			m1    = mfcba  - mfaba;
+			m0    = m2		+ mfbba;
+			mfaba = m0;
+			mfbba = m1 -		m0 * vvx;
+			mfcba = m2 - two*	m1 * vvx + vx2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfaca  + mfcca;
+			m1    = mfcca  - mfaca;
+			m0    = m2		+ mfbca;
+			mfaca = m0;
+			m0   += c1o3 * oMdrho;
+			mfbca = m1 -		m0 * vvx;
+			mfcca = m2 - two*	m1 * vvx + vx2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfaab	+ mfcab;
+			m1    = mfcab	- mfaab;
+			m0    = m2		+ mfbab;
+			mfaab = m0;
+			mfbab = m1 -		m0 * vvx;
+			mfcab = m2 - two*	m1 * vvx + vx2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfabb  + mfcbb;
+			m1    = mfcbb  - mfabb;
+			m0    = m2		+ mfbbb;
+			mfabb = m0;
+			mfbbb = m1 -		m0 * vvx;
+			mfcbb = m2 - two*	m1 * vvx + vx2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfacb  + mfccb;
+			m1    = mfccb  - mfacb;
+			m0    = m2		+ mfbcb;
+			mfacb = m0;
+			mfbcb = m1 -		m0 * vvx;
+			mfccb = m2 - two*	m1 * vvx + vx2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfaac	+ mfcac;
+			m1    = mfcac	- mfaac;
+			m0    = m2		+ mfbac;
+			mfaac = m0;
+			m0   += c1o3 * oMdrho;
+			mfbac = m1 -		m0 * vvx;
+			mfcac = m2 - two*	m1 * vvx + vx2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfabc  + mfcbc;
+			m1    = mfcbc  - mfabc;
+			m0    = m2		+ mfbbc;
+			mfabc = m0;
+			mfbbc = m1 -		m0 * vvx;
+			mfcbc = m2 - two*	m1 * vvx + vx2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			m2    = mfacc  + mfccc;
+			m1    = mfccc  - mfacc;
+			m0    = m2		+ mfbcc;
+			mfacc = m0;
+			m0   += c1o9 * oMdrho;
+			mfbcc = m1 -		m0 * vvx;
+			mfccc = m2 - two*	m1 * vvx + vx2 * m0;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+
+
+			////////////////////////////////////////////////////////////////////////////////////
+			// Cumulants
+			////////////////////////////////////////////////////////////////////////////////////
+			doubflo OxxPyyPzz = one;	//set the bulk viscosity one is high / two is very low and zero is (too) high
+
+			////////////////////////////////////////////////////////////
+			//3.
+			//////////////////////////////
+			doubflo OxyyPxzz  = one;
+			doubflo OxyyMxzz  = one;
+			doubflo Oxyz      = one;
+			////////////////////////////////////////////////////////////
+			//4.
+			//////////////////////////////
+			doubflo O4        = one;
+			////////////////////////////////////////////////////////////
+			//5.
+			//////////////////////////////
+			doubflo O5        = one;
+			////////////////////////////////////////////////////////////
+			//6.
+			//////////////////////////////
+			doubflo O6        = one;
+			////////////////////////////////////////////////////////////
+
+
+			//central moments to cumulants
+			//4.
+			doubflo CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + two * mfbba * mfbab) / rho;
+			doubflo CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + two * mfbba * mfabb) / rho;
+			doubflo CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + two * mfbab * mfabb) / rho;
+				  	 		
+			doubflo CUMcca = mfcca - (((mfcaa * mfaca + two * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho  - c1o9*(drho/rho));
+			doubflo CUMcac = mfcac - (((mfcaa * mfaac + two * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho  - c1o9*(drho/rho));
+			doubflo CUMacc = mfacc - (((mfaac * mfaca + two * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho  - c1o9*(drho/rho));
+
+			//5.
+			doubflo CUMbcc = mfbcc - ((mfaac * mfbca + mfaca * mfbac + four * mfabb * mfbbb + two * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac) ) / rho ;
+			doubflo CUMcbc = mfcbc - ((mfaac * mfcba + mfcaa * mfabc + four * mfbab * mfbbb + two * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc) ) / rho ;
+			doubflo CUMccb = mfccb - ((mfcaa * mfacb + mfaca * mfcab + four * mfbba * mfbbb + two * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab) ) / rho ;
+			
+			//6.
+
+			doubflo CUMccc = mfccc + ((-four *  mfbbb * mfbbb  
+							-           (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+							-    four * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+							-     two * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) / rho
+							+(   four * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+							+     two * (mfcaa * mfaca * mfaac)
+							+ sixteen *  mfbba * mfbab * mfabb) / (rho * rho)
+							-    c1o3 * (mfacc + mfcac + mfcca) /rho 
+							-    c1o9 * (mfcaa + mfaca + mfaac) /rho 
+							+(    two * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) 
+							+           (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) / (rho * rho) * c2o3 
+							+ c1o27*((drho * drho - drho)/(rho*rho)));
+							//+ c1o27*(one -three/rho +two/(rho*rho)));
+
+
+
+
+			//2.
+			// linear combinations
+			doubflo mxxPyyPzz = mfcaa + mfaca + mfaac;
+			doubflo mxxMyy    = mfcaa - mfaca;
+			doubflo mxxMzz	  = mfcaa - mfaac;
+			
+			//////////////////////////////////////////////////////////////////////////
+// 			doubflo magicBulk=(CUMacc+CUMcac+CUMcca)*(one/OxxPyyPzz-c1o2)*c3o2*8.;
+
+			//////////////////////////////////////////////////////////////////////////
+			//limiter-Scheise Teil 1
+			//doubflo oxxyy,oxxzz,oxy,oxz,oyz;
+			//doubflo smag=0.001;
+			//oxxyy    = omega+(one-omega)*abs(mxxMyy)/(abs(mxxMyy)+smag);
+			//oxxzz    = omega+(one-omega)*abs(mxxMzz)/(abs(mxxMzz)+smag);
+			//oxy      = omega+(one-omega)*abs(mfbba)/(abs(mfbba)+smag);
+			//oxz      = omega+(one-omega)*abs(mfbab)/(abs(mfbab)+smag);
+			//oyz      = omega+(one-omega)*abs(mfabb)/(abs(mfabb)+smag);
+
+			////////////////////////////////////////////////////////////////////////////
+			////Teil 1b
+			//doubflo constante = 1000.0;
+			//doubflo nuEddi = constante * abs(mxxPyyPzz);
+			//doubflo omegaLimit = one / (one / omega + three * nuEddi);
+
+			//{
+			//	doubflo dxux = c1o2 * (-omegaLimit) *(mxxMyy + mxxMzz) +  OxxPyyPzz * (mfaaa - mxxPyyPzz);
+			//	doubflo dyuy = dxux + omegaLimit * c3o2 * mxxMyy;
+			//	doubflo dzuz = dxux + omegaLimit * c3o2 * mxxMzz;
+
+				////relax
+				//mxxPyyPzz += OxxPyyPzz*(mfaaa  - mxxPyyPzz)- three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
+				//mxxMyy    += omegaLimit * (-mxxMyy) - three * (one + c1o2 * (-omegaLimit)) * (vx2 * dxux + vy2 * dyuy);
+				//mxxMzz    += omegaLimit * (-mxxMzz) - three * (one + c1o2 * (-omegaLimit)) * (vx2 * dxux + vz2 * dzuz);
+
+			//}
+			//mfabb     += omegaLimit * (-mfabb);
+			//mfbab     += omegaLimit * (-mfbab);
+			//mfbba     += omegaLimit * (-mfbba);
+			////////////////////////////////////////////////////////////////////////////
+
+ 			///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ 			//incl. correction		(hat noch nicht so gut funktioniert...Optimierungsbedarf??)
+ 			{
+ 				doubflo dxux = c1o2 * (-omega) *(mxxMyy + mxxMzz) + c1o2 *  OxxPyyPzz * (mfaaa - mxxPyyPzz);
+ 				doubflo dyuy = dxux + omega * c3o2 * mxxMyy;
+ 				doubflo dzuz = dxux + omega * c3o2 * mxxMzz;
+ 
+ 				////relax original
+ 				//mxxPyyPzz += OxxPyyPzz*(mfaaa  - mxxPyyPzz)- three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);//-magicBulk*OxxPyyPzz;
+ 				//mxxMyy    += omega * (-mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
+ 				//mxxMzz    += omega * (-mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
+				//relax porous media Mp
+				mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz + rho * ( vx2 + vy2 + vz2) * (one / porosity - one)) - three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);//-magicBulk*OxxPyyPzz;
+			    //relax porous media Ms
+				mxxMyy += omega * (rho * (vx2 - vy2) * (one / porosity - one) - mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
+				mxxMzz += omega * (rho * (vx2 - vz2) * (one / porosity - one) - mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
+
+ 				//////////////////////////////////////////////////////////////////////////
+ 				//limiter-Scheise Teil 2
+ 				//mxxMyy    += oxxyy * (-mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux + vy2 * dyuy);
+ 				//mxxMzz    += oxxzz * (-mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux + vz2 * dzuz);
+ 				//////////////////////////////////////////////////////////////////////////
+ 
+ 			}
+ 			///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ 			/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ 			////no correction
+ 			//mxxPyyPzz += OxxPyyPzz*(mfaaa-mxxPyyPzz);//-magicBulk*OxxPyyPzz;
+ 			//mxxMyy    += -(-omega) * (-mxxMyy);
+ 			//mxxMzz    += -(-omega) * (-mxxMzz);
+ 			/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			////original
+			//mfabb     += omega * (-mfabb);
+			//mfbab     += omega * (-mfbab);
+			//mfbba     += omega * (-mfbba);
+			/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			//porous media M11
+			mfabb += omega * (rho * vvy * vvz * (one / porosity - one) - mfabb);
+			mfbab += omega * (rho * vvx * vvz * (one / porosity - one) - mfbab);
+			mfbba += omega * (rho * vvx * vvy * (one / porosity - one) - mfbba);
+
+			//////////////////////////////////////////////////////////////////////////
+			//limiter-Scheise Teil 3
+			//mfabb     += oyz * (-mfabb);
+			//mfbab     += oxz * (-mfbab);
+			//mfbba     += oxy * (-mfbba);
+			//////////////////////////////////////////////////////////////////////////
+
+			// linear combinations back
+			mfcaa = c1o3 * (       mxxMyy +      mxxMzz + mxxPyyPzz);
+			mfaca = c1o3 * (-two*  mxxMyy +      mxxMzz + mxxPyyPzz);
+			mfaac = c1o3 * (       mxxMyy - two* mxxMzz + mxxPyyPzz);
+
+			//3.
+			// linear combinations
+
+			doubflo mxxyPyzz = mfcba + mfabc;
+			doubflo mxxyMyzz = mfcba - mfabc;
+
+			doubflo mxxzPyyz = mfcab + mfacb;
+			doubflo mxxzMyyz = mfcab - mfacb;
+
+			doubflo mxyyPxzz = mfbca + mfbac;
+			doubflo mxyyMxzz = mfbca - mfbac;
+
+			//relax
+			//////////////////////////////////////////////////////////////////////////
+			//das ist der limiter
+ 			wadjust    = Oxyz+(one-Oxyz)*abs(mfbbb)/(abs(mfbbb)+qudricLimitD);
+ 			mfbbb     += wadjust * (-mfbbb);
+ 			wadjust    = OxyyPxzz+(one-OxyyPxzz)*abs(mxxyPyzz)/(abs(mxxyPyzz)+qudricLimitP);
+ 			mxxyPyzz  += wadjust * (-mxxyPyzz);
+ 			wadjust    = OxyyMxzz+(one-OxyyMxzz)*abs(mxxyMyzz)/(abs(mxxyMyzz)+qudricLimitM);
+ 			mxxyMyzz  += wadjust * (-mxxyMyzz);
+ 			wadjust    = OxyyPxzz+(one-OxyyPxzz)*abs(mxxzPyyz)/(abs(mxxzPyyz)+qudricLimitP);
+ 			mxxzPyyz  += wadjust * (-mxxzPyyz);
+ 			wadjust    = OxyyMxzz+(one-OxyyMxzz)*abs(mxxzMyyz)/(abs(mxxzMyyz)+qudricLimitM);
+ 			mxxzMyyz  += wadjust * (-mxxzMyyz);
+ 			wadjust    = OxyyPxzz+(one-OxyyPxzz)*abs(mxyyPxzz)/(abs(mxyyPxzz)+qudricLimitP);
+ 			mxyyPxzz  += wadjust * (-mxyyPxzz);
+ 			wadjust    = OxyyMxzz+(one-OxyyMxzz)*abs(mxyyMxzz)/(abs(mxyyMxzz)+qudricLimitM);
+ 			mxyyMxzz  += wadjust * (-mxyyMxzz);
+			//////////////////////////////////////////////////////////////////////////
+			//ohne limiter
+			//mfbbb     += OxyyMxzz * (-mfbbb);
+			//mxxyPyzz  += OxyyPxzz * (-mxxyPyzz);
+			//mxxyMyzz  += OxyyMxzz * (-mxxyMyzz);
+			//mxxzPyyz  += OxyyPxzz * (-mxxzPyyz);
+			//mxxzMyyz  += OxyyMxzz * (-mxxzMyyz);
+			//mxyyPxzz  += OxyyPxzz * (-mxyyPxzz);
+			//mxyyMxzz  += OxyyMxzz * (-mxyyMxzz);
+			//////////////////////////////////////////////////////////////////////////
+
+			//// linear combinations back
+
+			mfcba = ( mxxyMyzz + mxxyPyzz) * c1o2;
+			mfabc = (-mxxyMyzz + mxxyPyzz) * c1o2;
+			mfcab = ( mxxzMyyz + mxxzPyyz) * c1o2;
+			mfacb = (-mxxzMyyz + mxxzPyyz) * c1o2;
+			mfbca = ( mxyyMxzz + mxyyPxzz) * c1o2;
+			mfbac = (-mxyyMxzz + mxyyPxzz) * c1o2;
+
+			//4.
+			//////////////////////////////////////////////////////////////////////////
+			//mit limiter
+ 		//	wadjust    = O4+(one-O4)*abs(CUMacc)/(abs(CUMacc)+qudricLimit);
+			//CUMacc    += wadjust * (-CUMacc);
+ 		//	wadjust    = O4+(one-O4)*abs(CUMcac)/(abs(CUMcac)+qudricLimit);
+			//CUMcac    += wadjust * (-CUMcac); 
+ 		//	wadjust    = O4+(one-O4)*abs(CUMcca)/(abs(CUMcca)+qudricLimit);
+			//CUMcca    += wadjust * (-CUMcca); 
+
+ 		//	wadjust    = O4+(one-O4)*abs(CUMbbc)/(abs(CUMbbc)+qudricLimit);
+			//CUMbbc    += wadjust * (-CUMbbc); 
+ 		//	wadjust    = O4+(one-O4)*abs(CUMbcb)/(abs(CUMbcb)+qudricLimit);
+			//CUMbcb    += wadjust * (-CUMbcb); 
+ 		//	wadjust    = O4+(one-O4)*abs(CUMcbb)/(abs(CUMcbb)+qudricLimit);
+			//CUMcbb    += wadjust * (-CUMcbb); 
+			//////////////////////////////////////////////////////////////////////////
+			//ohne limiter
+			CUMacc += O4 * (-CUMacc); 
+			CUMcac += O4 * (-CUMcac); 
+			CUMcca += O4 * (-CUMcca); 
+
+			CUMbbc += O4 * (-CUMbbc); 
+			CUMbcb += O4 * (-CUMbcb); 
+			CUMcbb += O4 * (-CUMcbb); 
+			//////////////////////////////////////////////////////////////////////////
+			
+					
+			//5.
+			CUMbcc += O5 * (-CUMbcc);
+			CUMcbc += O5 * (-CUMcbc);
+			CUMccb += O5 * (-CUMccb);
+
+			//6.
+			CUMccc += O6 * (-CUMccc);
+			
+
+
+			//back cumulants to central moments
+			//4.
+			mfcbb = CUMcbb + ((mfcaa + c1o3) * mfabb + two * mfbba * mfbab) / rho;
+			mfbcb = CUMbcb + ((mfaca + c1o3) * mfbab + two * mfbba * mfabb) / rho;
+			mfbbc = CUMbbc + ((mfaac + c1o3) * mfbba + two * mfbab * mfabb) / rho;
+						   
+			mfcca = CUMcca + (((mfcaa * mfaca + two * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho  - c1o9*(drho/rho));
+			mfcac = CUMcac + (((mfcaa * mfaac + two * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho  - c1o9*(drho/rho));
+			mfacc = CUMacc + (((mfaac * mfaca + two * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho  - c1o9*(drho/rho));
+
+			//5.
+			mfbcc = CUMbcc + ((mfaac * mfbca + mfaca * mfbac + four * mfabb * mfbbb + two * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac) ) / rho ;
+			mfcbc = CUMcbc + ((mfaac * mfcba + mfcaa * mfabc + four * mfbab * mfbbb + two * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc) ) / rho ;
+			mfccb = CUMccb + ((mfcaa * mfacb + mfaca * mfcab + four * mfbba * mfbbb + two * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab) ) / rho ;
+			
+			//6.
+
+			mfccc = CUMccc - ((-four *  mfbbb * mfbbb  
+							-           (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+							-    four * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+							-     two * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) / rho
+							+(   four * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+							+     two * (mfcaa * mfaca * mfaac)
+							+ sixteen *  mfbba * mfbab * mfabb) / (rho * rho)
+							-    c1o3 * (mfacc + mfcac + mfcca) /rho 
+							-    c1o9 * (mfcaa + mfaca + mfaac) /rho 
+							+(    two * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) 
+							+           (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) / (rho * rho) * c2o3 
+							+ c1o27*((drho * drho - drho)/(rho*rho)));
+							//+ c1o27*(one -three/rho +two/(rho*rho)));
+			////////////////////////////////////////////////////////////////////////////////////
+
+			////////////////////////////////////////////////////////////////////////////////////
+			//the force be with you
+			mfbaa = -mfbaa;
+			mfaba = -mfaba;
+			mfaab = -mfaab;
+			////////////////////////////////////////////////////////////////////////////////////
+
+
+			////////////////////////////////////////////////////////////////////////////////////
+			//back
+			////////////////////////////////////////////////////////////////////////////////////
+			//mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9   Konditionieren
+			////////////////////////////////////////////////////////////////////////////////////
+			// Z - Dir
+			m0 =  mfaac * c1o2 +      mfaab * (vvz - c1o2) + (mfaaa + one* oMdrho) * (     vz2 - vvz) * c1o2; 
+			m1 = -mfaac        - two* mfaab *  vvz         +  mfaaa                * (one- vz2)              - one* oMdrho * vz2; 
+			m2 =  mfaac * c1o2 +      mfaab * (vvz + c1o2) + (mfaaa + one* oMdrho) * (     vz2 + vvz) * c1o2;
+			mfaaa = m0;
+			mfaab = m1;
+			mfaac = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfabc * c1o2 +      mfabb * (vvz - c1o2) + mfaba * (     vz2 - vvz) * c1o2; 
+			m1 = -mfabc        - two* mfabb *  vvz         + mfaba * (one- vz2); 
+			m2 =  mfabc * c1o2 +      mfabb * (vvz + c1o2) + mfaba * (     vz2 + vvz) * c1o2;
+			mfaba = m0;
+			mfabb = m1;
+			mfabc = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfacc * c1o2 +      mfacb * (vvz - c1o2) + (mfaca + c1o3 * oMdrho) * (     vz2 - vvz) * c1o2; 
+			m1 = -mfacc        - two* mfacb *  vvz         +  mfaca                  * (one- vz2)              - c1o3 * oMdrho * vz2; 
+			m2 =  mfacc * c1o2 +      mfacb * (vvz + c1o2) + (mfaca + c1o3 * oMdrho) * (     vz2 + vvz) * c1o2;
+			mfaca = m0;
+			mfacb = m1;
+			mfacc = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfbac * c1o2 +      mfbab * (vvz - c1o2) + mfbaa * (     vz2 - vvz) * c1o2; 
+			m1 = -mfbac        - two* mfbab *  vvz         + mfbaa * (one- vz2); 
+			m2 =  mfbac * c1o2 +      mfbab * (vvz + c1o2) + mfbaa * (     vz2 + vvz) * c1o2;
+			mfbaa = m0;
+			mfbab = m1;
+			mfbac = m2;
+			/////////b//////////////////////////////////////////////////////////////////////////
+			m0 =  mfbbc * c1o2 +      mfbbb * (vvz - c1o2) + mfbba * (     vz2 - vvz) * c1o2; 
+			m1 = -mfbbc        - two* mfbbb *  vvz         + mfbba * (one- vz2); 
+			m2 =  mfbbc * c1o2 +      mfbbb * (vvz + c1o2) + mfbba * (     vz2 + vvz) * c1o2;
+			mfbba = m0;
+			mfbbb = m1;
+			mfbbc = m2;
+			/////////b//////////////////////////////////////////////////////////////////////////
+			m0 =  mfbcc * c1o2 +      mfbcb * (vvz - c1o2) + mfbca * (     vz2 - vvz) * c1o2; 
+			m1 = -mfbcc        - two* mfbcb *  vvz         + mfbca * (one- vz2); 
+			m2 =  mfbcc * c1o2 +      mfbcb * (vvz + c1o2) + mfbca * (     vz2 + vvz) * c1o2;
+			mfbca = m0;
+			mfbcb = m1;
+			mfbcc = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfcac * c1o2 +      mfcab * (vvz - c1o2) + (mfcaa + c1o3 * oMdrho) * (     vz2 - vvz) * c1o2; 
+			m1 = -mfcac        - two* mfcab *  vvz         +  mfcaa                  * (one- vz2)              - c1o3 * oMdrho * vz2; 
+			m2 =  mfcac * c1o2 +      mfcab * (vvz + c1o2) + (mfcaa + c1o3 * oMdrho) * (     vz2 + vvz) * c1o2;
+			mfcaa = m0;
+			mfcab = m1;
+			mfcac = m2;
+			/////////c//////////////////////////////////////////////////////////////////////////
+			m0 =  mfcbc * c1o2 +      mfcbb * (vvz - c1o2) + mfcba * (     vz2 - vvz) * c1o2; 
+			m1 = -mfcbc        - two* mfcbb *  vvz         + mfcba * (one- vz2); 
+			m2 =  mfcbc * c1o2 +      mfcbb * (vvz + c1o2) + mfcba * (     vz2 + vvz) * c1o2;
+			mfcba = m0;
+			mfcbb = m1;
+			mfcbc = m2;
+			/////////c//////////////////////////////////////////////////////////////////////////
+			m0 =  mfccc * c1o2 +      mfccb * (vvz - c1o2) + (mfcca + c1o9 * oMdrho) * (     vz2 - vvz) * c1o2; 
+			m1 = -mfccc        - two* mfccb *  vvz         +  mfcca                  * (one- vz2)              - c1o9 * oMdrho * vz2; 
+			m2 =  mfccc * c1o2 +      mfccb * (vvz + c1o2) + (mfcca + c1o9 * oMdrho) * (     vz2 + vvz) * c1o2;
+			mfcca = m0;
+			mfccb = m1;
+			mfccc = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			//mit 1/6, 2/3, 1/6, 0, 0, 0, 1/18, 2/9, 1/18   Konditionieren
+			////////////////////////////////////////////////////////////////////////////////////
+			// Y - Dir
+			m0 =  mfaca * c1o2 +      mfaba * (vvy - c1o2) + (mfaaa + c1o6 * oMdrho) * (     vy2 - vvy) * c1o2; 
+			m1 = -mfaca        - two* mfaba *  vvy         +  mfaaa                  * (one- vy2)              - c1o6 * oMdrho * vy2; 
+			m2 =  mfaca * c1o2 +      mfaba * (vvy + c1o2) + (mfaaa + c1o6 * oMdrho) * (     vy2 + vvy) * c1o2;
+			mfaaa = m0;
+			mfaba = m1;
+			mfaca = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfacb * c1o2 +      mfabb * (vvy - c1o2) + (mfaab + c2o3 * oMdrho) * (     vy2 - vvy) * c1o2; 
+			m1 = -mfacb        - two* mfabb *  vvy         +  mfaab                  * (one- vy2)              - c2o3 * oMdrho * vy2; 
+			m2 =  mfacb * c1o2 +      mfabb * (vvy + c1o2) + (mfaab + c2o3 * oMdrho) * (     vy2 + vvy) * c1o2;
+			mfaab = m0;
+			mfabb = m1;
+			mfacb = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfacc * c1o2 +      mfabc * (vvy - c1o2) + (mfaac + c1o6 * oMdrho) * (     vy2 - vvy) * c1o2; 
+			m1 = -mfacc        - two* mfabc *  vvy         +  mfaac                  * (one- vy2)              - c1o6 * oMdrho * vy2; 
+			m2 =  mfacc * c1o2 +      mfabc * (vvy + c1o2) + (mfaac + c1o6 * oMdrho) * (     vy2 + vvy) * c1o2;
+			mfaac = m0;
+			mfabc = m1;
+			mfacc = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfbca * c1o2 +      mfbba * (vvy - c1o2) + mfbaa * (     vy2 - vvy) * c1o2; 
+			m1 = -mfbca        - two* mfbba *  vvy         + mfbaa * (one- vy2); 
+			m2 =  mfbca * c1o2 +      mfbba * (vvy + c1o2) + mfbaa * (     vy2 + vvy) * c1o2;
+			mfbaa = m0;
+			mfbba = m1;
+			mfbca = m2;
+			/////////b//////////////////////////////////////////////////////////////////////////
+			m0 =  mfbcb * c1o2 +      mfbbb * (vvy - c1o2) + mfbab * (     vy2 - vvy) * c1o2; 
+			m1 = -mfbcb        - two* mfbbb *  vvy         + mfbab * (one- vy2); 
+			m2 =  mfbcb * c1o2 +      mfbbb * (vvy + c1o2) + mfbab * (     vy2 + vvy) * c1o2;
+			mfbab = m0;
+			mfbbb = m1;
+			mfbcb = m2;
+			/////////b//////////////////////////////////////////////////////////////////////////
+			m0 =  mfbcc * c1o2 +      mfbbc * (vvy - c1o2) + mfbac * (     vy2 - vvy) * c1o2; 
+			m1 = -mfbcc        - two* mfbbc *  vvy         + mfbac * (one- vy2); 
+			m2 =  mfbcc * c1o2 +      mfbbc * (vvy + c1o2) + mfbac * (     vy2 + vvy) * c1o2;
+			mfbac = m0;
+			mfbbc = m1;
+			mfbcc = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfcca * c1o2 +      mfcba * (vvy - c1o2) + (mfcaa + c1o18 * oMdrho) * (     vy2 - vvy) * c1o2; 
+			m1 = -mfcca        - two* mfcba *  vvy         +  mfcaa                   * (one- vy2)              - c1o18 * oMdrho * vy2; 
+			m2 =  mfcca * c1o2 +      mfcba * (vvy + c1o2) + (mfcaa + c1o18 * oMdrho) * (     vy2 + vvy) * c1o2;
+			mfcaa = m0;
+			mfcba = m1;
+			mfcca = m2;
+			/////////c//////////////////////////////////////////////////////////////////////////
+			m0 =  mfccb * c1o2 +      mfcbb * (vvy - c1o2) + (mfcab + c2o9 * oMdrho) * (     vy2 - vvy) * c1o2; 
+			m1 = -mfccb        - two* mfcbb *  vvy         +  mfcab                  * (one- vy2)              - c2o9 * oMdrho * vy2; 
+			m2 =  mfccb * c1o2 +      mfcbb * (vvy + c1o2) + (mfcab + c2o9 * oMdrho) * (     vy2 + vvy) * c1o2;
+			mfcab = m0;
+			mfcbb = m1;
+			mfccb = m2;
+			/////////c//////////////////////////////////////////////////////////////////////////
+			m0 =  mfccc * c1o2 +      mfcbc * (vvy - c1o2) + (mfcac + c1o18 * oMdrho) * (     vy2 - vvy) * c1o2; 
+			m1 = -mfccc        - two* mfcbc *  vvy         +  mfcac                   * (one- vy2)              - c1o18 * oMdrho * vy2; 
+			m2 =  mfccc * c1o2 +      mfcbc * (vvy + c1o2) + (mfcac + c1o18 * oMdrho) * (     vy2 + vvy) * c1o2;
+			mfcac = m0;
+			mfcbc = m1;
+			mfccc = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			//mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
+			////////////////////////////////////////////////////////////////////////////////////
+			// X - Dir
+			m0 =  mfcaa * c1o2 +      mfbaa * (vvx - c1o2) + (mfaaa + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
+			m1 = -mfcaa        - two* mfbaa *  vvx         +  mfaaa                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
+			m2 =  mfcaa * c1o2 +      mfbaa * (vvx + c1o2) + (mfaaa + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
+			mfaaa = m0;
+			mfbaa = m1;
+			mfcaa = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfcba * c1o2 +      mfbba * (vvx - c1o2) + (mfaba + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+			m1 = -mfcba        - two* mfbba *  vvx         +  mfaba                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
+			m2 =  mfcba * c1o2 +      mfbba * (vvx + c1o2) + (mfaba + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+			mfaba = m0;
+			mfbba = m1;
+			mfcba = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfcca * c1o2 +      mfbca * (vvx - c1o2) + (mfaca + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
+			m1 = -mfcca        - two* mfbca *  vvx         +  mfaca                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
+			m2 =  mfcca * c1o2 +      mfbca * (vvx + c1o2) + (mfaca + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
+			mfaca = m0;
+			mfbca = m1;
+			mfcca = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfcab * c1o2 +      mfbab * (vvx - c1o2) + (mfaab + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+			m1 = -mfcab        - two* mfbab *  vvx         +  mfaab                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
+			m2 =  mfcab * c1o2 +      mfbab * (vvx + c1o2) + (mfaab + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+			mfaab = m0;
+			mfbab = m1;
+			mfcab = m2;
+			///////////b////////////////////////////////////////////////////////////////////////
+			m0 =  mfcbb * c1o2 +      mfbbb * (vvx - c1o2) + (mfabb + c4o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+			m1 = -mfcbb        - two* mfbbb *  vvx         +  mfabb                  * (one- vx2)              - c4o9 * oMdrho * vx2; 
+			m2 =  mfcbb * c1o2 +      mfbbb * (vvx + c1o2) + (mfabb + c4o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+			mfabb = m0;
+			mfbbb = m1;
+			mfcbb = m2;
+			///////////b////////////////////////////////////////////////////////////////////////
+			m0 =  mfccb * c1o2 +      mfbcb * (vvx - c1o2) + (mfacb + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+			m1 = -mfccb        - two* mfbcb *  vvx         +  mfacb                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
+			m2 =  mfccb * c1o2 +      mfbcb * (vvx + c1o2) + (mfacb + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+			mfacb = m0;
+			mfbcb = m1;
+			mfccb = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+			////////////////////////////////////////////////////////////////////////////////////
+			m0 =  mfcac * c1o2 +      mfbac * (vvx - c1o2) + (mfaac + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
+			m1 = -mfcac        - two* mfbac *  vvx         +  mfaac                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
+			m2 =  mfcac * c1o2 +      mfbac * (vvx + c1o2) + (mfaac + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
+			mfaac = m0;
+			mfbac = m1;
+			mfcac = m2;
+			///////////c////////////////////////////////////////////////////////////////////////
+			m0 =  mfcbc * c1o2 +      mfbbc * (vvx - c1o2) + (mfabc + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+			m1 = -mfcbc        - two* mfbbc *  vvx         +  mfabc                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
+			m2 =  mfcbc * c1o2 +      mfbbc * (vvx + c1o2) + (mfabc + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+			mfabc = m0;
+			mfbbc = m1;
+			mfcbc = m2;
+			///////////c////////////////////////////////////////////////////////////////////////
+			m0 =  mfccc * c1o2 +      mfbcc * (vvx - c1o2) + (mfacc + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
+			m1 = -mfccc        - two* mfbcc *  vvx         +  mfacc                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
+			m2 =  mfccc * c1o2 +      mfbcc * (vvx + c1o2) + (mfacc + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
+			mfacc = m0;
+			mfbcc = m1;
+			mfccc = m2;
+			////////////////////////////////////////////////////////////////////////////////////
+
+			////////////////////////////////////////////////////////////////////////////////////
+			(D.f[ dirE   ])[k   ] = mfabb;                                                                     
+			(D.f[ dirW   ])[kw  ] = mfcbb;                                                                   
+			(D.f[ dirN   ])[k   ] = mfbab;
+			(D.f[ dirS   ])[ks  ] = mfbcb;
+			(D.f[ dirT   ])[k   ] = mfbba;
+			(D.f[ dirB   ])[kb  ] = mfbbc;
+			(D.f[ dirNE  ])[k   ] = mfaab;
+			(D.f[ dirSW  ])[ksw ] = mfccb;
+			(D.f[ dirSE  ])[ks  ] = mfacb;
+			(D.f[ dirNW  ])[kw  ] = mfcab;
+			(D.f[ dirTE  ])[k   ] = mfaba;
+			(D.f[ dirBW  ])[kbw ] = mfcbc;
+			(D.f[ dirBE  ])[kb  ] = mfabc;
+			(D.f[ dirTW  ])[kw  ] = mfcba;
+			(D.f[ dirTN  ])[k   ] = mfbaa;
+			(D.f[ dirBS  ])[kbs ] = mfbcc;
+			(D.f[ dirBN  ])[kb  ] = mfbac;
+			(D.f[ dirTS  ])[ks  ] = mfbca;
+			(D.f[ dirZERO])[k   ] = mfbbb;
+			(D.f[ dirTNE ])[k   ] = mfaaa;
+			(D.f[ dirTSE ])[ks  ] = mfaca;
+			(D.f[ dirBNE ])[kb  ] = mfaac;
+			(D.f[ dirBSE ])[kbs ] = mfacc;
+			(D.f[ dirTNW ])[kw  ] = mfcaa;
+			(D.f[ dirTSW ])[ksw ] = mfcca;
+			(D.f[ dirBNW ])[kbw ] = mfcac;
+			(D.f[ dirBSW ])[kbsw] = mfccc;
+			////////////////////////////////////////////////////////////////////////////////////
+		}                                                                                                                    
+	}
+}
+////////////////////////////////////////////////////////////////////////////////
+
diff --git a/src/VirtualFluids_GPU/GPU/VelocityBCs27.cu b/src/VirtualFluids_GPU/GPU/VelocityBCs27.cu
index a0f5b84218ae97468421586c78c242ec0853f8ea..6d88f4b5a4b55823e177254720d05ef27ee2bbe7 100644
--- a/src/VirtualFluids_GPU/GPU/VelocityBCs27.cu
+++ b/src/VirtualFluids_GPU/GPU/VelocityBCs27.cu
@@ -4964,7 +4964,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c2over27* (drho/*+three*( vx1        )*/+c9over2*( vx1        )*( vx1        ) * (one + drho)-cu_sq); 
-         (D.f[dirW])[kw]=(one-q)/(one+q)*(f_E-f_W+(f_E+f_W-two*feq*om1)/(one-om1))*c1o2+(q*(f_E+f_W)-six*c2over27*( VeloX     ))/(one+q);// - c2over27 * drho;
+         (D.f[dirW])[kw]=(one-q)/(one+q)*(f_E-f_W+(f_E+f_W-two*feq*om1)/(one-om1))*c1o2+(q*(f_E+f_W)-six*c2over27*( VeloX     ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
          //(D.f[dirW])[kw]=zero;
       }
 
@@ -4972,7 +4972,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c2over27* (drho/*+three*(-vx1        )*/+c9over2*(-vx1        )*(-vx1        ) * (one + drho)-cu_sq); 
-         (D.f[dirE])[ke]=(one-q)/(one+q)*(f_W-f_E+(f_W+f_E-two*feq*om1)/(one-om1))*c1o2+(q*(f_W+f_E)-six*c2over27*(-VeloX     ))/(one+q);// - c2over27 * drho;
+         (D.f[dirE])[ke]=(one-q)/(one+q)*(f_W-f_E+(f_W+f_E-two*feq*om1)/(one-om1))*c1o2+(q*(f_W+f_E)-six*c2over27*(-VeloX     ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
          //(D.f[dirE])[ke]=zero;
       }
 
@@ -4980,7 +4980,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c2over27* (drho/*+three*(    vx2     )*/+c9over2*(     vx2    )*(     vx2    ) * (one + drho)-cu_sq); 
-         (D.f[dirS])[ks]=(one-q)/(one+q)*(f_N-f_S+(f_N+f_S-two*feq*om1)/(one-om1))*c1o2+(q*(f_N+f_S)-six*c2over27*( VeloY     ))/(one+q);// - c2over27 * drho;
+         (D.f[dirS])[ks]=(one-q)/(one+q)*(f_N-f_S+(f_N+f_S-two*feq*om1)/(one-om1))*c1o2+(q*(f_N+f_S)-six*c2over27*( VeloY     ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
          //(D.f[dirS])[ks]=zero;
       }
 
@@ -4988,7 +4988,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c2over27* (drho/*+three*(   -vx2     )*/+c9over2*(    -vx2    )*(    -vx2    ) * (one + drho)-cu_sq); 
-         (D.f[dirN])[kn]=(one-q)/(one+q)*(f_S-f_N+(f_S+f_N-two*feq*om1)/(one-om1))*c1o2+(q*(f_S+f_N)-six*c2over27*(-VeloY     ))/(one+q);// - c2over27 * drho;
+         (D.f[dirN])[kn]=(one-q)/(one+q)*(f_S-f_N+(f_S+f_N-two*feq*om1)/(one-om1))*c1o2+(q*(f_S+f_N)-six*c2over27*(-VeloY     ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
          //(D.f[dirN])[kn]=zero;
       }
 
@@ -4996,7 +4996,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c2over27* (drho/*+three*(         vx3)*/+c9over2*(         vx3)*(         vx3) * (one + drho)-cu_sq); 
-         (D.f[dirB])[kb]=(one-q)/(one+q)*(f_T-f_B+(f_T+f_B-two*feq*om1)/(one-om1))*c1o2+(q*(f_T+f_B)-six*c2over27*( VeloZ     ))/(one+q);// - c2over27 * drho;
+         (D.f[dirB])[kb]=(one-q)/(one+q)*(f_T-f_B+(f_T+f_B-two*feq*om1)/(one-om1))*c1o2+(q*(f_T+f_B)-six*c2over27*( VeloZ     )/* * (one + drho)*/)/(one+q);// - c2over27 * drho;
          //(D.f[dirB])[kb]=one;
       }
 
@@ -5004,7 +5004,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c2over27* (drho/*+three*(        -vx3)*/+c9over2*(        -vx3)*(        -vx3) * (one + drho)-cu_sq); 
-         (D.f[dirT])[kt]=(one-q)/(one+q)*(f_B-f_T+(f_B+f_T-two*feq*om1)/(one-om1))*c1o2+(q*(f_B+f_T)-six*c2over27*(-VeloZ     ))/(one+q);// - c2over27 * drho;
+         (D.f[dirT])[kt]=(one-q)/(one+q)*(f_B-f_T+(f_B+f_T-two*feq*om1)/(one-om1))*c1o2+(q*(f_B+f_T)-six*c2over27*(-VeloZ     ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
          //(D.f[dirT])[kt]=zero;
       }
 
@@ -5012,7 +5012,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*( vx1+vx2    )*/+c9over2*( vx1+vx2    )*( vx1+vx2    ) * (one + drho)-cu_sq); 
-         (D.f[dirSW])[ksw]=(one-q)/(one+q)*(f_NE-f_SW+(f_NE+f_SW-two*feq*om1)/(one-om1))*c1o2+(q*(f_NE+f_SW)-six*c1over54*(VeloX+VeloY))/(one+q);// - c1over54 * drho;
+         (D.f[dirSW])[ksw]=(one-q)/(one+q)*(f_NE-f_SW+(f_NE+f_SW-two*feq*om1)/(one-om1))*c1o2+(q*(f_NE+f_SW)-six*c1over54*(VeloX+VeloY) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirSW])[ksw]=zero;
       }
 
@@ -5020,7 +5020,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*(-vx1-vx2    )*/+c9over2*(-vx1-vx2    )*(-vx1-vx2    ) * (one + drho)-cu_sq); 
-         (D.f[dirNE])[kne]=(one-q)/(one+q)*(f_SW-f_NE+(f_SW+f_NE-two*feq*om1)/(one-om1))*c1o2+(q*(f_SW+f_NE)-six*c1over54*(-VeloX-VeloY))/(one+q);// - c1over54 * drho;
+         (D.f[dirNE])[kne]=(one-q)/(one+q)*(f_SW-f_NE+(f_SW+f_NE-two*feq*om1)/(one-om1))*c1o2+(q*(f_SW+f_NE)-six*c1over54*(-VeloX-VeloY) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirNE])[kne]=zero;
       }
 
@@ -5028,7 +5028,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*( vx1-vx2    )*/+c9over2*( vx1-vx2    )*( vx1-vx2    ) * (one + drho)-cu_sq); 
-         (D.f[dirNW])[knw]=(one-q)/(one+q)*(f_SE-f_NW+(f_SE+f_NW-two*feq*om1)/(one-om1))*c1o2+(q*(f_SE+f_NW)-six*c1over54*( VeloX-VeloY))/(one+q);// - c1over54 * drho;
+         (D.f[dirNW])[knw]=(one-q)/(one+q)*(f_SE-f_NW+(f_SE+f_NW-two*feq*om1)/(one-om1))*c1o2+(q*(f_SE+f_NW)-six*c1over54*( VeloX-VeloY) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirNW])[knw]=zero;
       }
 
@@ -5036,7 +5036,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*(-vx1+vx2    )*/+c9over2*(-vx1+vx2    )*(-vx1+vx2    ) * (one + drho)-cu_sq); 
-         (D.f[dirSE])[kse]=(one-q)/(one+q)*(f_NW-f_SE+(f_NW+f_SE-two*feq*om1)/(one-om1))*c1o2+(q*(f_NW+f_SE)-six*c1over54*(-VeloX+VeloY))/(one+q);// - c1over54 * drho;
+         (D.f[dirSE])[kse]=(one-q)/(one+q)*(f_NW-f_SE+(f_NW+f_SE-two*feq*om1)/(one-om1))*c1o2+(q*(f_NW+f_SE)-six*c1over54*(-VeloX+VeloY) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirSE])[kse]=zero;
       }
 
@@ -5044,7 +5044,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*( vx1    +vx3)*/+c9over2*( vx1    +vx3)*( vx1    +vx3) * (one + drho)-cu_sq); 
-         (D.f[dirBW])[kbw]=(one-q)/(one+q)*(f_TE-f_BW+(f_TE+f_BW-two*feq*om1)/(one-om1))*c1o2+(q*(f_TE+f_BW)-six*c1over54*( VeloX+VeloZ))/(one+q);// - c1over54 * drho;
+         (D.f[dirBW])[kbw]=(one-q)/(one+q)*(f_TE-f_BW+(f_TE+f_BW-two*feq*om1)/(one-om1))*c1o2+(q*(f_TE+f_BW)-six*c1over54*( VeloX+VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirBW])[kbw]=zero;
       }
 
@@ -5052,7 +5052,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*(-vx1    -vx3)*/+c9over2*(-vx1    -vx3)*(-vx1    -vx3) * (one + drho)-cu_sq); 
-         (D.f[dirTE])[kte]=(one-q)/(one+q)*(f_BW-f_TE+(f_BW+f_TE-two*feq*om1)/(one-om1))*c1o2+(q*(f_BW+f_TE)-six*c1over54*(-VeloX-VeloZ))/(one+q);// - c1over54 * drho;
+         (D.f[dirTE])[kte]=(one-q)/(one+q)*(f_BW-f_TE+(f_BW+f_TE-two*feq*om1)/(one-om1))*c1o2+(q*(f_BW+f_TE)-six*c1over54*(-VeloX-VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirTE])[kte]=zero;
       }
 
@@ -5060,7 +5060,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*( vx1    -vx3)*/+c9over2*( vx1    -vx3)*( vx1    -vx3) * (one + drho)-cu_sq); 
-         (D.f[dirTW])[ktw]=(one-q)/(one+q)*(f_BE-f_TW+(f_BE+f_TW-two*feq*om1)/(one-om1))*c1o2+(q*(f_BE+f_TW)-six*c1over54*( VeloX-VeloZ))/(one+q);// - c1over54 * drho;
+         (D.f[dirTW])[ktw]=(one-q)/(one+q)*(f_BE-f_TW+(f_BE+f_TW-two*feq*om1)/(one-om1))*c1o2+(q*(f_BE+f_TW)-six*c1over54*( VeloX-VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirTW])[ktw]=zero;
       }
 
@@ -5068,7 +5068,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*(-vx1    +vx3)*/+c9over2*(-vx1    +vx3)*(-vx1    +vx3) * (one + drho)-cu_sq); 
-         (D.f[dirBE])[kbe]=(one-q)/(one+q)*(f_TW-f_BE+(f_TW+f_BE-two*feq*om1)/(one-om1))*c1o2+(q*(f_TW+f_BE)-six*c1over54*(-VeloX+VeloZ))/(one+q);// - c1over54 * drho;
+         (D.f[dirBE])[kbe]=(one-q)/(one+q)*(f_TW-f_BE+(f_TW+f_BE-two*feq*om1)/(one-om1))*c1o2+(q*(f_TW+f_BE)-six*c1over54*(-VeloX+VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirBE])[kbe]=zero;
       }
 
@@ -5076,7 +5076,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*(     vx2+vx3)*/+c9over2*(     vx2+vx3)*(     vx2+vx3) * (one + drho)-cu_sq); 
-         (D.f[dirBS])[kbs]=(one-q)/(one+q)*(f_TN-f_BS+(f_TN+f_BS-two*feq*om1)/(one-om1))*c1o2+(q*(f_TN+f_BS)-six*c1over54*( VeloY+VeloZ))/(one+q);// - c1over54 * drho;
+         (D.f[dirBS])[kbs]=(one-q)/(one+q)*(f_TN-f_BS+(f_TN+f_BS-two*feq*om1)/(one-om1))*c1o2+(q*(f_TN+f_BS)-six*c1over54*( VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirBS])[kbs]=zero;
       }
 
@@ -5084,7 +5084,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*(    -vx2-vx3)*/+c9over2*(    -vx2-vx3)*(    -vx2-vx3) * (one + drho)-cu_sq); 
-         (D.f[dirTN])[ktn]=(one-q)/(one+q)*(f_BS-f_TN+(f_BS+f_TN-two*feq*om1)/(one-om1))*c1o2+(q*(f_BS+f_TN)-six*c1over54*( -VeloY-VeloZ))/(one+q);// - c1over54 * drho;
+         (D.f[dirTN])[ktn]=(one-q)/(one+q)*(f_BS-f_TN+(f_BS+f_TN-two*feq*om1)/(one-om1))*c1o2+(q*(f_BS+f_TN)-six*c1over54*( -VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirTN])[ktn]=zero;
       }
 
@@ -5092,7 +5092,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*(     vx2-vx3)*/+c9over2*(     vx2-vx3)*(     vx2-vx3) * (one + drho)-cu_sq); 
-         (D.f[dirTS])[kts]=(one-q)/(one+q)*(f_BN-f_TS+(f_BN+f_TS-two*feq*om1)/(one-om1))*c1o2+(q*(f_BN+f_TS)-six*c1over54*( VeloY-VeloZ))/(one+q);// - c1over54 * drho;
+         (D.f[dirTS])[kts]=(one-q)/(one+q)*(f_BN-f_TS+(f_BN+f_TS-two*feq*om1)/(one-om1))*c1o2+(q*(f_BN+f_TS)-six*c1over54*( VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirTS])[kts]=zero;
       }
 
@@ -5100,7 +5100,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over54* (drho/*+three*(    -vx2+vx3)*/+c9over2*(    -vx2+vx3)*(    -vx2+vx3) * (one + drho)-cu_sq); 
-         (D.f[dirBN])[kbn]=(one-q)/(one+q)*(f_TS-f_BN+(f_TS+f_BN-two*feq*om1)/(one-om1))*c1o2+(q*(f_TS+f_BN)-six*c1over54*( -VeloY+VeloZ))/(one+q);// - c1over54 * drho;
+         (D.f[dirBN])[kbn]=(one-q)/(one+q)*(f_TS-f_BN+(f_TS+f_BN-two*feq*om1)/(one-om1))*c1o2+(q*(f_TS+f_BN)-six*c1over54*( -VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
          //(D.f[dirBN])[kbn]=zero;
       }
 
@@ -5108,7 +5108,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over216*(drho/*+three*( vx1+vx2+vx3)*/+c9over2*( vx1+vx2+vx3)*( vx1+vx2+vx3) * (one + drho)-cu_sq); 
-         (D.f[dirBSW])[kbsw]=(one-q)/(one+q)*(f_TNE-f_BSW+(f_TNE+f_BSW-two*feq*om1)/(one-om1))*c1o2+(q*(f_TNE+f_BSW)-six*c1over216*( VeloX+VeloY+VeloZ))/(one+q);// - c1over216 * drho;
+         (D.f[dirBSW])[kbsw]=(one-q)/(one+q)*(f_TNE-f_BSW+(f_TNE+f_BSW-two*feq*om1)/(one-om1))*c1o2+(q*(f_TNE+f_BSW)-six*c1over216*( VeloX+VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
          //(D.f[dirBSW])[kbsw]=zero;
       }
 
@@ -5116,7 +5116,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over216*(drho/*+three*(-vx1-vx2-vx3)*/+c9over2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3) * (one + drho)-cu_sq); 
-         (D.f[dirTNE])[ktne]=(one-q)/(one+q)*(f_BSW-f_TNE+(f_BSW+f_TNE-two*feq*om1)/(one-om1))*c1o2+(q*(f_BSW+f_TNE)-six*c1over216*(-VeloX-VeloY-VeloZ))/(one+q);// - c1over216 * drho;
+         (D.f[dirTNE])[ktne]=(one-q)/(one+q)*(f_BSW-f_TNE+(f_BSW+f_TNE-two*feq*om1)/(one-om1))*c1o2+(q*(f_BSW+f_TNE)-six*c1over216*(-VeloX-VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
          //(D.f[dirTNE])[ktne]=zero;
       }
 
@@ -5124,7 +5124,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over216*(drho/*+three*( vx1+vx2-vx3)*/+c9over2*( vx1+vx2-vx3)*( vx1+vx2-vx3) * (one + drho)-cu_sq); 
-         (D.f[dirTSW])[ktsw]=(one-q)/(one+q)*(f_BNE-f_TSW+(f_BNE+f_TSW-two*feq*om1)/(one-om1))*c1o2+(q*(f_BNE+f_TSW)-six*c1over216*( VeloX+VeloY-VeloZ))/(one+q);// - c1over216 * drho;
+         (D.f[dirTSW])[ktsw]=(one-q)/(one+q)*(f_BNE-f_TSW+(f_BNE+f_TSW-two*feq*om1)/(one-om1))*c1o2+(q*(f_BNE+f_TSW)-six*c1over216*( VeloX+VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
          //(D.f[dirTSW])[ktsw]=zero;
       }
 
@@ -5132,7 +5132,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over216*(drho/*+three*(-vx1-vx2+vx3)*/+c9over2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3) * (one + drho)-cu_sq); 
-         (D.f[dirBNE])[kbne]=(one-q)/(one+q)*(f_TSW-f_BNE+(f_TSW+f_BNE-two*feq*om1)/(one-om1))*c1o2+(q*(f_TSW+f_BNE)-six*c1over216*(-VeloX-VeloY+VeloZ))/(one+q);// - c1over216 * drho;
+         (D.f[dirBNE])[kbne]=(one-q)/(one+q)*(f_TSW-f_BNE+(f_TSW+f_BNE-two*feq*om1)/(one-om1))*c1o2+(q*(f_TSW+f_BNE)-six*c1over216*(-VeloX-VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
          //(D.f[dirBNE])[kbne]=zero;
       }
 
@@ -5140,7 +5140,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over216*(drho/*+three*( vx1-vx2+vx3)*/+c9over2*( vx1-vx2+vx3)*( vx1-vx2+vx3) * (one + drho)-cu_sq); 
-         (D.f[dirBNW])[kbnw]=(one-q)/(one+q)*(f_TSE-f_BNW+(f_TSE+f_BNW-two*feq*om1)/(one-om1))*c1o2+(q*(f_TSE+f_BNW)-six*c1over216*( VeloX-VeloY+VeloZ))/(one+q);// - c1over216 * drho;
+         (D.f[dirBNW])[kbnw]=(one-q)/(one+q)*(f_TSE-f_BNW+(f_TSE+f_BNW-two*feq*om1)/(one-om1))*c1o2+(q*(f_TSE+f_BNW)-six*c1over216*( VeloX-VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
          //(D.f[dirBNW])[kbnw]=zero;
       }
 
@@ -5148,7 +5148,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over216*(drho/*+three*(-vx1+vx2-vx3)*/+c9over2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3) * (one + drho)-cu_sq); 
-         (D.f[dirTSE])[ktse]=(one-q)/(one+q)*(f_BNW-f_TSE+(f_BNW+f_TSE-two*feq*om1)/(one-om1))*c1o2+(q*(f_BNW+f_TSE)-six*c1over216*(-VeloX+VeloY-VeloZ))/(one+q);// - c1over216 * drho;
+         (D.f[dirTSE])[ktse]=(one-q)/(one+q)*(f_BNW-f_TSE+(f_BNW+f_TSE-two*feq*om1)/(one-om1))*c1o2+(q*(f_BNW+f_TSE)-six*c1over216*(-VeloX+VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
          //(D.f[dirTSE])[ktse]=zero;
       }
 
@@ -5156,7 +5156,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over216*(drho/*+three*( vx1-vx2-vx3)*/+c9over2*( vx1-vx2-vx3)*( vx1-vx2-vx3) * (one + drho)-cu_sq); 
-         (D.f[dirTNW])[ktnw]=(one-q)/(one+q)*(f_BSE-f_TNW+(f_BSE+f_TNW-two*feq*om1)/(one-om1))*c1o2+(q*(f_BSE+f_TNW)-six*c1over216*( VeloX-VeloY-VeloZ))/(one+q);// - c1over216 * drho;
+         (D.f[dirTNW])[ktnw]=(one-q)/(one+q)*(f_BSE-f_TNW+(f_BSE+f_TNW-two*feq*om1)/(one-om1))*c1o2+(q*(f_BSE+f_TNW)-six*c1over216*( VeloX-VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
          //(D.f[dirTNW])[ktnw]=zero;
       }
 
@@ -5164,7 +5164,7 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
       if (q>=zero && q<=one)
       {
          feq=c1over216*(drho/*+three*(-vx1+vx2+vx3)*/+c9over2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3) * (one + drho)-cu_sq); 
-         (D.f[dirBSE])[kbse]=(one-q)/(one+q)*(f_TNW-f_BSE+(f_TNW+f_BSE-two*feq*om1)/(one-om1))*c1o2+(q*(f_TNW+f_BSE)-six*c1over216*(-VeloX+VeloY+VeloZ))/(one+q);// - c1over216 * drho;
+         (D.f[dirBSE])[kbse]=(one-q)/(one+q)*(f_TNW-f_BSE+(f_TNW+f_BSE-two*feq*om1)/(one-om1))*c1o2+(q*(f_TNW+f_BSE)-six*c1over216*(-VeloX+VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
          //(D.f[dirBSE])[kbse]=zero;
       }
    }
diff --git a/src/VirtualFluids_GPU/Interface_OpenFOAM/Interface.cpp b/src/VirtualFluids_GPU/Interface_OpenFOAM/Interface.cpp
index edb304f9bf197571660b36212852bc084bcaedff..a9009cb73e178384915f70a556c88b4d428ac48e 100644
--- a/src/VirtualFluids_GPU/Interface_OpenFOAM/Interface.cpp
+++ b/src/VirtualFluids_GPU/Interface_OpenFOAM/Interface.cpp
@@ -981,7 +981,7 @@ void Interface::allocArrays_BoundaryValues(Parameter* para) {
 				//para->getParH(i)->Qinflow.Vx[m] = para->getParH(i)->Qinflow.Vx[m] / para->getVelocityRatio();
 				//para->getParH(i)->Qinflow.Vy[m] = para->getParH(i)->Qinflow.Vy[m] / para->getVelocityRatio();
 				//para->getParH(i)->Qinflow.Vz[m] = para->getParH(i)->Qinflow.Vz[m] / para->getVelocityRatio();
-				para->getParH(i)->Qinflow.Vx[m] = 0.0; //para->getVelocity();//0.035;
+				para->getParH(i)->Qinflow.Vx[m] = para->getVelocity();//0.035;
 				para->getParH(i)->Qinflow.Vy[m] = 0.0;//para->getVelocity();//0.0;
 				para->getParH(i)->Qinflow.Vz[m] = 0.0;
 				//if (para->getParH(i)->Qinflow.Vz[m] > 0)
diff --git a/src/VirtualFluids_GPU/LBM/LB.h b/src/VirtualFluids_GPU/LBM/LB.h
index fa89a202962d3418c1b5b2a2900d5d29175d1e4f..e454f8468752e8a420cf1835614cf1f5f3714a9f 100644
--- a/src/VirtualFluids_GPU/LBM/LB.h
+++ b/src/VirtualFluids_GPU/LBM/LB.h
@@ -6,6 +6,13 @@
 #define GEO_VELO         2
 #define GEO_PRESS        4
 
+//////////////////////////
+//porous media
+#define GEO_PM_0		 5
+#define GEO_PM_1		 6
+#define GEO_PM_2		 7
+//////////////////////////
+
 #define GEO_SOLID       15
 #define GEO_VOID        16
 
diff --git a/src/VirtualFluids_GPU/LBM/Simulation.cpp b/src/VirtualFluids_GPU/LBM/Simulation.cpp
index d19e4e8e4837e1d0b230a22e0e83efa2a8dbd62b..1569eda2598183c6a8d41c101841c17117813cdf 100644
--- a/src/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -37,6 +37,7 @@
 #include "Calculation/Calc2ndMoments.h"
 #include "Calculation/CalcMedian.h"
 #include "Calculation/ForceCalculations.h"
+#include "Calculation/PorousMedia.h"
 //////////////////////////////////////////////////////////////////////////
 //CUDA
 #include <cuda_runtime.h>
@@ -544,79 +545,101 @@ void Simulation::run()
 		 //getLastCudaError("KernelCasSPKum27 execution failed");
 		 //////////////////////////////////////////////////////////////////////////
 
-		//////////////////////////////////////////////////////////////////////////
-		//Wale 
-		if (para->getUseWale())
-		{
-			KernelWaleCumOneCompSP27(para->getParD(0)->numberofthreads,
-									 para->getParD(0)->omega,			
-									 para->getParD(0)->geoSP, 
-									 para->getParD(0)->neighborX_SP, 
-									 para->getParD(0)->neighborY_SP, 
-									 para->getParD(0)->neighborZ_SP,
-									 para->getParD(0)->neighborWSB_SP,
-								     para->getParD(0)->vx_SP,        
-								     para->getParD(0)->vy_SP,        
-								     para->getParD(0)->vz_SP,        
-									 para->getParD(0)->d0SP.f[0],
-									 para->getParD(0)->turbViscosity,
-									 para->getParD(0)->size_Mat_SP,
-									 para->getParD(0)->size_Array_SP,
-									 0,
-									 para->getForcesDev(),
-									 para->getParD(0)->evenOrOdd); 
-			getLastCudaError("KernelWaleCumOneCompSP27 execution failed");
-
-			//KernelWaleCumAA2016CompSP27(para->getParD(0)->numberofthreads,
-			//							para->getParD(0)->omega,			
-			//							para->getParD(0)->geoSP, 
-			//							para->getParD(0)->neighborX_SP, 
-			//							para->getParD(0)->neighborY_SP, 
-			//							para->getParD(0)->neighborZ_SP,
-			//							para->getParD(0)->neighborWSB_SP,
-			//							para->getParD(0)->vx_SP,        
-			//							para->getParD(0)->vy_SP,        
-			//							para->getParD(0)->vz_SP,        
-			//							para->getParD(0)->d0SP.f[0],
-			//							para->getParD(0)->turbViscosity,
-			//							para->getParD(0)->size_Mat_SP,
-			//							para->getParD(0)->size_Array_SP,
-			//							0,
-			//							para->getForcesDev(),
-			//							para->getParD(0)->evenOrOdd); 
-			//getLastCudaError("KernelWaleCumAA2016CompSP27 execution failed");
+		////////////////////////////////////////////////////////////////////////////
+		////Wale 
+		//if (para->getUseWale())
+		//{
+		//	KernelWaleCumOneCompSP27(para->getParD(0)->numberofthreads,
+		//							 para->getParD(0)->omega,			
+		//							 para->getParD(0)->geoSP, 
+		//							 para->getParD(0)->neighborX_SP, 
+		//							 para->getParD(0)->neighborY_SP, 
+		//							 para->getParD(0)->neighborZ_SP,
+		//							 para->getParD(0)->neighborWSB_SP,
+		//						     para->getParD(0)->vx_SP,        
+		//						     para->getParD(0)->vy_SP,        
+		//						     para->getParD(0)->vz_SP,        
+		//							 para->getParD(0)->d0SP.f[0],
+		//							 para->getParD(0)->turbViscosity,
+		//							 para->getParD(0)->size_Mat_SP,
+		//							 para->getParD(0)->size_Array_SP,
+		//							 0,
+		//							 para->getForcesDev(),
+		//							 para->getParD(0)->evenOrOdd); 
+		//	getLastCudaError("KernelWaleCumOneCompSP27 execution failed");
+
+		//	//KernelWaleCumAA2016CompSP27(para->getParD(0)->numberofthreads,
+		//	//							para->getParD(0)->omega,			
+		//	//							para->getParD(0)->geoSP, 
+		//	//							para->getParD(0)->neighborX_SP, 
+		//	//							para->getParD(0)->neighborY_SP, 
+		//	//							para->getParD(0)->neighborZ_SP,
+		//	//							para->getParD(0)->neighborWSB_SP,
+		//	//							para->getParD(0)->vx_SP,        
+		//	//							para->getParD(0)->vy_SP,        
+		//	//							para->getParD(0)->vz_SP,        
+		//	//							para->getParD(0)->d0SP.f[0],
+		//	//							para->getParD(0)->turbViscosity,
+		//	//							para->getParD(0)->size_Mat_SP,
+		//	//							para->getParD(0)->size_Array_SP,
+		//	//							0,
+		//	//							para->getForcesDev(),
+		//	//							para->getParD(0)->evenOrOdd); 
+		//	//getLastCudaError("KernelWaleCumAA2016CompSP27 execution failed");
+
+		//} 
+		//else
+		//{
+		//	KernelKumNewCompSP27(para->getParD(0)->numberofthreads,       
+		//						 para->getParD(0)->omega,			
+		//						 para->getParD(0)->geoSP, 
+		//						 para->getParD(0)->neighborX_SP, 
+		//						 para->getParD(0)->neighborY_SP, 
+		//						 para->getParD(0)->neighborZ_SP,
+		//						 para->getParD(0)->d0SP.f[0],    
+		//						 para->getParD(0)->size_Mat_SP,
+		//						 para->getParD(0)->size_Array_SP,
+		//						 0,
+		//						 para->getForcesDev(),
+		//						 para->getParD(0)->evenOrOdd); 
+		//	getLastCudaError("KernelCasSPKum27 execution failed");
+
+		//	//KernelKumAA2016CompSP27(para->getParD(0)->numberofthreads,       
+		//	//					 para->getParD(0)->omega, 
+		//	//					 para->getParD(0)->geoSP, 
+		//	//					 para->getParD(0)->neighborX_SP, 
+		//	//					 para->getParD(0)->neighborY_SP, 
+		//	//					 para->getParD(0)->neighborZ_SP,
+		//	//					 para->getParD(0)->d0SP.f[0],    
+		//	//					 para->getParD(0)->size_Mat_SP,
+		//	//					 0,
+		//	//					 para->getForcesDev(),
+		//	//					 para->getParD(0)->evenOrOdd); 
+		//	//getLastCudaError("KernelKumAA2016CompSP27 execution failed");
+
+		//}
+		////////////////////////////////////////////////////////////////////////////
 
-		} 
-		else
-		{
-			KernelKumNewCompSP27(para->getParD(0)->numberofthreads,       
-								 para->getParD(0)->omega,			
-								 para->getParD(0)->geoSP, 
-								 para->getParD(0)->neighborX_SP, 
-								 para->getParD(0)->neighborY_SP, 
-								 para->getParD(0)->neighborZ_SP,
-								 para->getParD(0)->d0SP.f[0],    
-								 para->getParD(0)->size_Mat_SP,
-								 para->getParD(0)->size_Array_SP,
-								 0,
-								 para->getForcesDev(),
-								 para->getParD(0)->evenOrOdd); 
-			getLastCudaError("KernelCasSPKum27 execution failed");
-
-			//KernelKumAA2016CompSP27(para->getParD(0)->numberofthreads,       
-			//					 para->getParD(0)->omega, 
-			//					 para->getParD(0)->geoSP, 
-			//					 para->getParD(0)->neighborX_SP, 
-			//					 para->getParD(0)->neighborY_SP, 
-			//					 para->getParD(0)->neighborZ_SP,
-			//					 para->getParD(0)->d0SP.f[0],    
-			//					 para->getParD(0)->size_Mat_SP,
-			//					 0,
-			//					 para->getForcesDev(),
-			//					 para->getParD(0)->evenOrOdd); 
-			//getLastCudaError("KernelKumAA2016CompSP27 execution failed");
 
-		}
+		//////////////////////////////////////////////////////////////////////////
+		//porous media
+		KernelPMCumOneCompSP27( para->getParD(0)->numberofthreads,
+								para->getParD(0)->omega,			
+								para->getParD(0)->geoSP, 
+								para->getParD(0)->neighborX_SP, 
+								para->getParD(0)->neighborY_SP, 
+								para->getParD(0)->neighborZ_SP,
+								para->getParD(0)->d0SP.f[0],    
+								para->getParD(0)->size_Mat_SP,
+								0,
+								para->getForcesDev(),
+								0.5,
+								134.4,
+								0.0068287,
+								GEO_PM_1,
+								para->getParD(0)->evenOrOdd); 
+		getLastCudaError("KernelPMCumOneCompSP27 execution failed");
+
 
 
 		 //////////////////////////////////////////////////////////////////////////
@@ -1045,13 +1068,13 @@ void Simulation::run()
 		    //            para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
 		    //            para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
 		    //  getLastCudaError("QVelDev27 execution failed");
-		     // QVelDevComp27(para->getParD(0)->numberofthreads, para->getParD(0)->nx,           para->getParD(0)->ny,
-							//para->getParD(0)->Qinflow.Vx,      para->getParD(0)->Qinflow.Vy,   para->getParD(0)->Qinflow.Vz,
-							//para->getParD(0)->d0SP.f[0],       para->getParD(0)->Qinflow.k,    para->getParD(0)->Qinflow.q27[0], 
-							//para->getParD(0)->kInflowQ,        para->getParD(0)->kInflowQ,     para->getParD(0)->omega,
-							//para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
-							//para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
-		     // getLastCudaError("QVelDevComp27 execution failed");
+		      QVelDevComp27(para->getParD(0)->numberofthreads, para->getParD(0)->nx,           para->getParD(0)->ny,
+							para->getParD(0)->Qinflow.Vx,      para->getParD(0)->Qinflow.Vy,   para->getParD(0)->Qinflow.Vz,
+							para->getParD(0)->d0SP.f[0],       para->getParD(0)->Qinflow.k,    para->getParD(0)->Qinflow.q27[0], 
+							para->getParD(0)->kInflowQ,        para->getParD(0)->kInflowQ,     para->getParD(0)->omega,
+							para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
+							para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
+		      getLastCudaError("QVelDevComp27 execution failed");
 		     // QVelDevComp27(para->getParD(0)->numberofthreads, para->getParD(0)->nx,           para->getParD(0)->ny,
 							//para->getParD(0)->QGeom.Vx,        para->getParD(0)->QGeom.Vy,     para->getParD(0)->QGeom.Vz,
 							//para->getParD(0)->d0SP.f[0],       para->getParD(0)->QGeom.k,      para->getParD(0)->QGeom.q27[0], 
@@ -1264,15 +1287,15 @@ void Simulation::run()
 			//		para->getParD(0)->neighborX_SP,          para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
 			//		para->getParD(0)->size_Mat_SP,           para->getParD(0)->evenOrOdd);
 			//getLastCudaError("QDev27 execution failed");
-		   if (para->getParD(0)->kQ > 0)
-		   {
-			   QDevComp27(para->getParD(0)->numberofthreads,       para->getParD(0)->nx,           para->getParD(0)->ny,
-			   		      para->getParD(0)->d0SP.f[0],             para->getParD(0)->QWall.k,	   para->getParD(0)->QWall.q27[0], 
-			   		      para->getParD(0)->kQ,                    para->getParD(0)->kQ,           para->getParD(0)->omega,
-			   		      para->getParD(0)->neighborX_SP,          para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
-			   		      para->getParD(0)->size_Mat_SP,           para->getParD(0)->evenOrOdd);
-			   getLastCudaError("QDevComp27 (Wall) execution failed");
-		   }
+		   //if (para->getParD(0)->kQ > 0)
+		   //{
+			  // QDevComp27(para->getParD(0)->numberofthreads,       para->getParD(0)->nx,           para->getParD(0)->ny,
+			  // 		      para->getParD(0)->d0SP.f[0],             para->getParD(0)->QWall.k,	   para->getParD(0)->QWall.q27[0], 
+			  // 		      para->getParD(0)->kQ,                    para->getParD(0)->kQ,           para->getParD(0)->omega,
+			  // 		      para->getParD(0)->neighborX_SP,          para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
+			  // 		      para->getParD(0)->size_Mat_SP,           para->getParD(0)->evenOrOdd);
+			  // getLastCudaError("QDevComp27 (Wall) execution failed");
+		   //}
 			////////////////////////////////////////////////////////////////////////////
 			//Geom
       //      BBDev27( para->getParD(0)->numberofthreads,       para->getParD(0)->nx,           para->getParD(0)->ny,
@@ -1364,11 +1387,11 @@ void Simulation::run()
 
 		} 
 		  //////////////////////////////////////////////////////////////////////////////////
-		  //calculate the new forcing
-		  if (((t%10) == 0) && (t >= 10)/*((t%para->getTStartOut()) == 0) && (t >= para->getTStartOut())*/)
-		  {
-			  forceCalculator->calcPIDControllerForForce(para);
-		  }
+		  ////calculate the new forcing
+		  //if (((t%10) == 0) && (t >= 10)/*((t%para->getTStartOut()) == 0) && (t >= para->getTStartOut())*/)
+		  //{
+			 // forceCalculator->calcPIDControllerForForce(para);
+		  //}
 		  //////////////////////////////////////////////////////////////////////////////////
 		  
 		  // ////////////////////////////////////////////////////////////////////////
@@ -1397,19 +1420,19 @@ void Simulation::run()
 		 //if (para->getParD(0)->QPress.kQ > 0)
 		 //{
 			// // //////////////////////////////////////////////////////////////////////////////////
-			// //QPressNoRhoDev27(  para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
-			// //					para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
-			// //					para->getParD(0)->QPress.kN,       para->getParD(0)->QPress.kQ,    para->getParD(0)->omega,
-			// //					para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
-			// //					para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
-			// //getLastCudaError("QPressNoRhoDev27 execution failed");
-			// //QPressDevEQZ27(para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
-		 //	//				para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
-			//	//			para->getParD(0)->QPress.kN,       para->getParD(0)->kDistTestRE.f[0],       
-			//	//			para->getParD(0)->QPress.kQ,       para->getParD(0)->omega,
-			//	//			para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
-		 //	//				para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
-			// //getLastCudaError("QPressDevEQZ27 execution failed");
+			 //QPressNoRhoDev27(  para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
+			 //					para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
+			 //					para->getParD(0)->QPress.kN,       para->getParD(0)->QPress.kQ,    para->getParD(0)->omega,
+			 //					para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
+			 //					para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
+			 //getLastCudaError("QPressNoRhoDev27 execution failed");
+			 //QPressDevEQZ27(para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
+		 	//				para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
+				//			para->getParD(0)->QPress.kN,       para->getParD(0)->kDistTestRE.f[0],       
+				//			para->getParD(0)->QPress.kQ,       para->getParD(0)->omega,
+				//			para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
+		 	//				para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
+			 //getLastCudaError("QPressDevEQZ27 execution failed");
 
  		////	QInflowScaleByPressDev27(   para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
 			////							para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
@@ -1442,12 +1465,12 @@ void Simulation::run()
 		  //getLastCudaError("QPressDevIncompNEQ27 execution failed");
 		  //////////////////////////////////////////////////////////////////////////////////
 		  //press NEQ comp
-		  //QPressDevNEQ27( para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
-		  //				  para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
-		  //				  para->getParD(0)->QPress.kN,       para->getParD(0)->QPress.kQ,    para->getParD(0)->omega,
-		  //				  para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
-		  //				  para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
-		  //getLastCudaError("QPressDevNEQ27 execution failed");
+		  QPressDevNEQ27( para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
+		  				  para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
+		  				  para->getParD(0)->QPress.kN,       para->getParD(0)->QPress.kQ,    para->getParD(0)->omega,
+		  				  para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
+		  				  para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
+		  getLastCudaError("QPressDevNEQ27 execution failed");
 		  ////////////////////////////////////////////////////////////////////////////////
           //if (  myid == numprocs - 1)  
           //   PressSchlaffer27( para->getParD(0)->numberofthreads,  para->getParD(0)->Qoutflow.RhoBC,
@@ -2582,3 +2605,10 @@ void Simulation::run()
    }
    //////////////////////////////////////////////////////////////////////////
 }
+void Simulation::porousMedia()
+{
+	//Kondensator = porous media 0
+	double porosity = 0.5;
+	//pm0 = new PorousMedia(porosity, GEO_PM_0,)
+
+}
diff --git a/src/VirtualFluids_GPU/LBM/Simulation.h b/src/VirtualFluids_GPU/LBM/Simulation.h
index b4dacb32493c048cae2f76c2711cd2bde2e297be..7fc3539d9ecc60668087022db78fd45d4143d17e 100644
--- a/src/VirtualFluids_GPU/LBM/Simulation.h
+++ b/src/VirtualFluids_GPU/LBM/Simulation.h
@@ -6,6 +6,7 @@
 #include "Utilities/Buffer2D.hpp"
 //#include "Input/ConfigFile.h"
 #include "Calculation/ForceCalculations.h"
+#include "Calculation/PorousMedia.h"
 #include "Parameter/Parameter.h"
 #include "Restart/RestartPostprocessor.h"
 #include "Restart/RestartObject.h"
@@ -22,6 +23,7 @@ public:
 	void run();
 	void init(std::string &cstr);
 	void bulk();
+	void porousMedia();
 protected:
 
 	Buffer2D <doubflo> sbuf_t; 
@@ -48,6 +50,11 @@ protected:
 	//Forcing Calculation
 	ForceCalculations* forceCalculator;
 
+	//Porous Media
+	PorousMedia* pm0;
+	PorousMedia* pm1;
+	PorousMedia* pm2;
+
 	//KQ - Schlaff
 	unsigned int            kNQ, kSQ, kEQ, kWQ;
 	QforBoundaryConditions  QnH, QnD;
diff --git a/src/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/VirtualFluids_GPU/Parameter/Parameter.cpp
index 6adf1e9af282b6e96e752b8a8e19bffb3daad461..bca9b54d3eee82e28ddf25152d18f4b24aab1952 100644
--- a/src/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -1966,6 +1966,40 @@ void Parameter::cudaAllocRandomValues()
 	checkCudaErrors( cudaMalloc((void**)&(this->devState), (sizeof(curandState)*parD[getFine()]->plp.numberOfParticles) ));
 }
 //////////////////////////////////////////////////////////////////////////
+//porous media
+void Parameter::cudaAllocPorousMedia(PorousMedia* pm, int pmID, int lev)
+{
+	unsigned int mem_size_IDsPM = sizeof(unsigned int)*pm[pmID].getSizePM();
+	unsigned int *tmpIDHost, *tmpIDDevice;
+
+	//Host
+	checkCudaErrors(cudaMallocHost((void**) &(tmpIDHost), mem_size_IDsPM));
+
+	//Device
+	checkCudaErrors(cudaMalloc((void**) &(tmpIDDevice), mem_size_IDsPM));
+
+	//////////////////////////////////////////////////////////////////////////
+	pm[pmID].setHostNodeIDsPM(tmpIDHost);
+	pm[pmID].setDeviceNodeIDsPM(tmpIDDevice);
+	//////////////////////////////////////////////////////////////////////////
+	double tmp = (double)mem_size_IDsPM;
+	setMemsizeGPU(tmp, false);
+}
+void Parameter::cudaCopyPorousMedia(PorousMedia* pm, int pmID, int lev)
+{
+	unsigned int mem_size_IDsPM = sizeof(unsigned int)*pm[pmID].getSizePM();
+	unsigned int *tmpIDHost   = pm[pmID].getHostNodeIDsPM();
+	unsigned int *tmpIDDevice = pm[pmID].getDeviceNodeIDsPM();
+	//////////////////////////////////////////////////////////////////////////
+	checkCudaErrors(cudaMemcpy(tmpIDDevice, tmpIDHost, mem_size_IDsPM, cudaMemcpyHostToDevice));
+	//////////////////////////////////////////////////////////////////////////
+	pm[pmID].setDeviceNodeIDsPM(tmpIDDevice);
+}
+void Parameter::cudaFreePorousMedia(PorousMedia* pm, int pmID, int lev)
+{
+	checkCudaErrors(cudaFreeHost(pm[pmID].getHostNodeIDsPM()));
+}
+//////////////////////////////////////////////////////////////////////////
 //advection diffusion
 void Parameter::cudaAllocConc(int lev)
 {
diff --git a/src/VirtualFluids_GPU/Parameter/Parameter.h b/src/VirtualFluids_GPU/Parameter/Parameter.h
index f7360b383973589267b167f7105d3ae45ba2a39e..63c3f2d019a7e49124aa6dd9314f0a180395a722 100644
--- a/src/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/VirtualFluids_GPU/Parameter/Parameter.h
@@ -4,6 +4,7 @@
 #include <vector>
 #include "LBM/LB.h"
 #include "LBM/D3Q27.h"
+#include "Calculation/PorousMedia.h"
 //#include "Output/LogWriter.hpp"
 //#include "boost/serialization/serialization.hpp"
 //#include "boost/serialization/vector.hpp"
@@ -450,6 +451,14 @@ public:
 	//////////////////////////////////////////////////////////////////////////
 
 
+	//////////////////////////////////////////////////////////////////////////
+	//Porous Media
+	void cudaAllocPorousMedia(PorousMedia* pm, int pmID, int lev);
+	void cudaCopyPorousMedia(PorousMedia* pm, int pmID, int lev);
+	void cudaFreePorousMedia(PorousMedia* pm, int pmID, int lev);
+	//////////////////////////////////////////////////////////////////////////
+
+
 	//////////////////////////////////////////////////////////////////////////
 	//Temperature
 	void cudaAllocConc(int lev);