diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 07ff3f6e3f200a4a2cde0d6fc6499b0cd8abb23d..7d3669b49163626b9a05fb2cf845f07c871439e9 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -526,11 +526,13 @@ doxygen:
   script:
   - apk update && apk add doxygen
   - doxygen docs/Doxyfile
+  - mv docs/build/html/ public/
 
   artifacts:
     expire_in: 1 week
     paths:
-    - docs/build/
+    - docs/build/html/
+    - public
 
 
 ###############################################################################
diff --git a/apps/gpu/LBM/WTG_RUB/CMakeLists.txt b/apps/gpu/LBM/WTG_RUB/CMakeLists.txt
index 8bac54bf335a26aaecbf8cbb5a1474e67bd99b89..606987dfb093c9c93bbd25bf5ff68fdc81311e1b 100644
--- a/apps/gpu/LBM/WTG_RUB/CMakeLists.txt
+++ b/apps/gpu/LBM/WTG_RUB/CMakeLists.txt
@@ -3,3 +3,5 @@ PROJECT(WTG_RUB LANGUAGES CUDA CXX)
 vf_add_library(BUILDTYPE binary PRIVATE_LINK basics VirtualFluids_GPU GridGenerator MPI::MPI_CXX FILES WTG_RUB.cpp)
 
 set_source_files_properties(WTG_RUB.cpp PROPERTIES LANGUAGE CUDA)
+
+set_target_properties(WTG_RUB PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
diff --git a/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp b/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp
index 12bc850bc718952278b5f0e9e10470eea3ff18b3..2528cd00a7e3b016767158eae085af596d6ea956 100644
--- a/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp
+++ b/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp
@@ -28,6 +28,7 @@
 #include "GridGenerator/grid/GridBuilder/LevelGridBuilder.h"
 #include "GridGenerator/grid/GridBuilder/MultipleGridBuilder.h"
 #include "GridGenerator/grid/BoundaryConditions/Side.h"
+#include "GridGenerator/grid/BoundaryConditions/BoundaryCondition.h"
 #include "GridGenerator/grid/GridFactory.h"
 #include "GridGenerator/geometries/TriangularMesh/TriangularMesh.h"
 #include "GridGenerator/geometries/Conglomerate/Conglomerate.h"
@@ -67,26 +68,24 @@ LbmOrGks lbmOrGks = LBM;
 
 const real L  = 1.0;
 
-const real Re = 500.0;// 1000.0;
-
 const real velocity  = 1.0;
 
 int variant = 1;
 real rotationOfCity;
 real z_offset = 0.0; // only if baseplate is in use (currently not!! not important)
+int dataN;
+std::vector<real> dataZ;
+std::vector<real> dataVelocity;
+std::string simulationName("");
 
 // 1: original setup of Lennard Lux (6 level, 4.0 cm -> 1.25 mm)
 // 2: setup 1 of MSch               (4 level, 1.0 cm -> 1.25 mm)
 // 3: setup 2 of MSch               (5 level, 1.6 cm -> 1.0  mm)
 int setupDomain = 3;
 
-
-
 std::string path("D:/out/WTG_RUB"); //Mollok
 std::string inputPath("D:/out/WTG_RUB/input/");
 
-std::string simulationName("RUB");
-
 const uint timeStepStartOut = 0;
 const uint timeStepOut = 10000;
 const uint timeStepEnd = 100000;
@@ -96,6 +95,8 @@ const uint timeStepEnd = 100000;
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 void addFineGrids(SPtr<MultipleGridBuilder> gridBuilder, uint &maxLevel, real &rotationOfCity);
 
+void readVelocityProfile();
+
 std::string chooseVariation();
 
 void multipleLevel(const std::string& configPath)
@@ -149,6 +150,26 @@ void multipleLevel(const std::string& configPath)
     TriangularMesh *RubSTL      = TriangularMesh::make(inputPath + "stl/" + chooseVariation() + ".stl");
     // vector<real> originOfCityXY = { 600.0, y_max / 2, z_offset };
 
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // OPTIONS
+
+    int maxLevelAtInflow = maxLevel;
+    // there are max. 3 levels (0,1,2) at inflow-/x_min-Side if maxLevel >= 3; OTHERWISE: maxLevel
+    if (maxLevel >= 3)
+        maxLevelAtInflow = 2;
+
+    // MeasurePoints [MP01-15: lvl maxLevel],[MP16-41: lvl 1]; disable when reducing numberOfLevels --> dx might be too
+    // large if MP01-15 are used with low resolution dx, MPs might be placed in solid City-geometry
+    bool useMP                   = true;
+    bool measureVeloProfilesOnly = false;
+
+    // Two Components: true->DiffOn, false->DiffOff
+    bool diffOnOff = false;
+
+    // Resetting diff or flow field, e.g. after restart, do not reset diff/flow at start of measureRun ;-)
+    bool reset_diff = true;
+    bool reset_flow = true;
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     gridBuilder->addCoarseGrid(x_min, y_min, z_min, 
                                x_max, y_max, z_max, dx);
@@ -225,6 +246,67 @@ void multipleLevel(const std::string& configPath)
 
     gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
 
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    readVelocityProfile();
+    for (int level = 0; level <= maxLevelAtInflow; level++) {
+        auto inflowBC = std::dynamic_pointer_cast<VelocityBoundaryCondition>(
+            gridBuilder->getBoundaryCondition(SideType::MX, level));
+
+        // lambda function [capture list](parameter list){function body}
+        // for every z, loop is being run;
+        if (inflowBC) {
+            inflowBC->setVelocityProfile(
+                gridBuilder->getGrid(level), [&](real x, real y, real z, real &vx, real &vy, real &vz) {
+                    int i;
+                    for (i = 0; i < dataN; i++) {
+                        if ((z < dataZ[i]) || (z == dataZ[i]))
+                            break;
+                    } // remembers i
+
+                    // if z is below current data point --> interpolation between previous and current datapoint
+                    if (z < dataZ[i]) {
+                        vx = velocityLB * (dataVelocity[i] + (dataVelocity[i + 1] - dataVelocity[i]) /
+                                                                 (dataZ[i + 1] - dataZ[i]) * (z - dataZ[i]));
+                    } else if (z == dataZ[i]) {
+                        vx = velocityLB * dataVelocity[i];
+                    }
+
+                    // vx << std::endl; vx = velocityLB;
+                    vy = 0.0;
+                    vz = 0.0;
+                });
+        }
+    }
+
+    // Resetting Velocity Profile (returning to intial state; do not combine with restart)
+    if (reset_flow) {
+        para->setInitialCondition([&](real x, real y, real z, real &rho, real &vx, real &vy, real &vz) {
+            int i;
+            for (i = 0; i < dataN; i++) {
+                if ((z < dataZ[i]) || (z == dataZ[i]))
+                    break;
+            }
+            if (z < dataZ[i]) {
+                vx = velocityLB * 0.588 *
+                     (dataVelocity[i] +
+                      (dataVelocity[i + 1] - dataVelocity[i]) / (dataZ[i + 1] - dataZ[i]) * (z - dataZ[i]));
+            } else if (z == dataZ[i]) {
+                vx = velocityLB * 0.588 * dataVelocity[i];
+            }
+
+            // std::cout << "LINE : " << __LINE__ << "\tdataZ = " << dataZ[i] << "\tz = " << z << "\tvx = " << vx <<
+            // std::endl;
+
+            // vx = velocityLB;
+            vy  = 0.0;
+            vz  = 0.0;
+            rho = 0.0;
+        });
+    }
+
+
+
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     SPtr<CudaMemoryManager> cudaMemoryManager = CudaMemoryManager::make(para);
@@ -242,6 +324,35 @@ void multipleLevel(const std::string& configPath)
 
 }
 
+void readVelocityProfile()
+{
+    // reads velocityProfile.txt, containing values for relative velocity (u_h) in relation to height (z). also fills
+    // dataZ,dataVelocity vectors
+
+    std::ifstream inFile;
+    inFile.open(inputPath + "VeloProfile.txt");
+    //inFile.open(inputPath + "velocityProfile.txt");
+
+    if (inFile.fail()) {
+        std::cerr << "Error opening File" << std::endl;
+        exit(1);
+    }
+
+    int z;
+    real velocity;
+
+    // read
+    inFile >> dataN;
+    for (int k = 0; k < dataN; k++) {
+        inFile >> z;
+        dataZ.push_back(z + z_offset);
+
+        inFile >> velocity;
+        dataVelocity.push_back(velocity);
+    }
+    inFile.close();
+}
+
 void addFineGrids(SPtr<MultipleGridBuilder> gridBuilder, uint &maxLevel, real &rotationOfCity)
 {
     if (setupDomain == 1) {
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu b/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
index 547b20745d4a6ab7209bb253b5a37c4eda13602e..ce6e034d1c1eee57e062b736cfcea97e07306f3c 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
@@ -42,14 +42,18 @@ extern "C" __global__ void LBCalcMac27( real* vxD,
 
    const unsigned int k = nx*(ny*z + y) + x; // Zugriff auf arrays im device
 
+
+   if(k >= size_Mat)
+      return;
+
+   if(!vf::gpu::isValidFluidNode(geoD[k]))
+      return;
+
    rhoD[k] = c0o1;
    vxD[k]  = c0o1;
    vyD[k]  = c0o1;
    vzD[k]  = c0o1;
 
-   if(!vf::gpu::isValidFluidNode(k, size_Mat, geoD[k]))
-      return;
-
    vf::gpu::DistributionWrapper distr_wrapper(distributions, size_Mat, isEvenTimestep, k, neighborX, neighborY, neighborZ);
    const auto& distribution = distr_wrapper.distribution;
 
@@ -264,206 +268,29 @@ extern "C" __global__ void LBCalcMacCompSP27(real *vxD, real *vyD, real *vzD, re
                                              unsigned int *neighborZ, unsigned int size_Mat, real *distributions,
                                              bool isEvenTimestep)
 {
-    //const unsigned k = vf::gpu::getNodeIndex();
-
-    //pressD[k] = c0o1;
-    //rhoD[k]   = c0o1;
-    //vxD[k]    = c0o1;
-    //vyD[k]    = c0o1;
-    //vzD[k]    = c0o1;
+    const unsigned k = vf::gpu::getNodeIndex();
 
-    //if (!vf::gpu::isValidFluidNode(k, size_Mat, geoD[k]))
-    //    return;
+    if(k >= size_Mat)
+        return;
 
-    //vf::gpu::DistributionWrapper distr_wrapper(distributions, size_Mat, isEvenTimestep, k, neighborX, neighborY,
-    //                                           neighborZ);
-    //const auto &distribution = distr_wrapper.distribution;
+    if (!vf::gpu::isValidFluidNode(geoD[k]))
+        return;
 
-    //rhoD[k]   = vf::lbm::getDensity(distribution.f);
-    //vxD[k]    = vf::lbm::getCompressibleVelocityX1(distribution.f, rhoD[k]);
-    //vyD[k]    = vf::lbm::getCompressibleVelocityX2(distribution.f, rhoD[k]);
-    //vzD[k]    = vf::lbm::getCompressibleVelocityX3(distribution.f, rhoD[k]);
-    //pressD[k] = vf::lbm::getPressure(distribution.f, rhoD[k], vxD[k], vyD[k], vzD[k]); 
+    pressD[k] = c0o1;
+    rhoD[k]   = c0o1;
+    vxD[k]    = c0o1;
+    vyD[k]    = c0o1;
+    vzD[k]    = c0o1;
 
+    vf::gpu::DistributionWrapper distr_wrapper(distributions, size_Mat, isEvenTimestep, k, neighborX, neighborY,
+                                               neighborZ);
+    const auto &distribution = distr_wrapper.distribution;
 
-   // old stuff
-   Distributions27 D;
-    if (isEvenTimestep == true)
-   {
-      D.f[dirE   ] = &distributions[dirE   *size_Mat];
-      D.f[dirW   ] = &distributions[dirW   *size_Mat];
-      D.f[dirN   ] = &distributions[dirN   *size_Mat];
-      D.f[dirS   ] = &distributions[dirS   *size_Mat];
-      D.f[dirT   ] = &distributions[dirT   *size_Mat];
-      D.f[dirB   ] = &distributions[dirB   *size_Mat];
-      D.f[dirNE  ] = &distributions[dirNE  *size_Mat];
-      D.f[dirSW  ] = &distributions[dirSW  *size_Mat];
-      D.f[dirSE  ] = &distributions[dirSE  *size_Mat];
-      D.f[dirNW  ] = &distributions[dirNW  *size_Mat];
-      D.f[dirTE  ] = &distributions[dirTE  *size_Mat];
-      D.f[dirBW  ] = &distributions[dirBW  *size_Mat];
-      D.f[dirBE  ] = &distributions[dirBE  *size_Mat];
-      D.f[dirTW  ] = &distributions[dirTW  *size_Mat];
-      D.f[dirTN  ] = &distributions[dirTN  *size_Mat];
-      D.f[dirBS  ] = &distributions[dirBS  *size_Mat];
-      D.f[dirBN  ] = &distributions[dirBN  *size_Mat];
-      D.f[dirTS  ] = &distributions[dirTS  *size_Mat];
-      D.f[dirZERO] = &distributions[dirZERO*size_Mat];
-      D.f[dirTNE ] = &distributions[dirTNE *size_Mat];
-      D.f[dirTSW ] = &distributions[dirTSW *size_Mat];
-      D.f[dirTSE ] = &distributions[dirTSE *size_Mat];
-      D.f[dirTNW ] = &distributions[dirTNW *size_Mat];
-      D.f[dirBNE ] = &distributions[dirBNE *size_Mat];
-      D.f[dirBSW ] = &distributions[dirBSW *size_Mat];
-      D.f[dirBSE ] = &distributions[dirBSE *size_Mat];
-      D.f[dirBNW ] = &distributions[dirBNW *size_Mat];
-   } 
-   else
-   {
-      D.f[dirW   ] = &distributions[dirE   *size_Mat];
-      D.f[dirE   ] = &distributions[dirW   *size_Mat];
-      D.f[dirS   ] = &distributions[dirN   *size_Mat];
-      D.f[dirN   ] = &distributions[dirS   *size_Mat];
-      D.f[dirB   ] = &distributions[dirT   *size_Mat];
-      D.f[dirT   ] = &distributions[dirB   *size_Mat];
-      D.f[dirSW  ] = &distributions[dirNE  *size_Mat];
-      D.f[dirNE  ] = &distributions[dirSW  *size_Mat];
-      D.f[dirNW  ] = &distributions[dirSE  *size_Mat];
-      D.f[dirSE  ] = &distributions[dirNW  *size_Mat];
-      D.f[dirBW  ] = &distributions[dirTE  *size_Mat];
-      D.f[dirTE  ] = &distributions[dirBW  *size_Mat];
-      D.f[dirTW  ] = &distributions[dirBE  *size_Mat];
-      D.f[dirBE  ] = &distributions[dirTW  *size_Mat];
-      D.f[dirBS  ] = &distributions[dirTN  *size_Mat];
-      D.f[dirTN  ] = &distributions[dirBS  *size_Mat];
-      D.f[dirTS  ] = &distributions[dirBN  *size_Mat];
-      D.f[dirBN  ] = &distributions[dirTS  *size_Mat];
-      D.f[dirZERO] = &distributions[dirZERO*size_Mat];
-      D.f[dirTNE ] = &distributions[dirBSW *size_Mat];
-      D.f[dirTSW ] = &distributions[dirBNE *size_Mat];
-      D.f[dirTSE ] = &distributions[dirBNW *size_Mat];
-      D.f[dirTNW ] = &distributions[dirBSE *size_Mat];
-      D.f[dirBNE ] = &distributions[dirTSW *size_Mat];
-      D.f[dirBSW ] = &distributions[dirTNE *size_Mat];
-      D.f[dirBSE ] = &distributions[dirTNW *size_Mat];
-      D.f[dirBNW ] = &distributions[dirTSE *size_Mat];
-   }
-   ////////////////////////////////////////////////////////////////////////////////
-   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)
-   {
-      //////////////////////////////////////////////////////////////////////////
-      //index
-      unsigned int kzero= k;
-      unsigned int ke   = k;
-      unsigned int kw   = neighborX[k];
-      unsigned int kn   = k;
-      unsigned int ks   = neighborY[k];
-      unsigned int kt   = k;
-      unsigned int kb   = neighborZ[k];
-      unsigned int ksw  = neighborY[kw];
-      unsigned int kne  = k;
-      unsigned int kse  = ks;
-      unsigned int knw  = kw;
-      unsigned int kbw  = neighborZ[kw];
-      unsigned int kte  = k;
-      unsigned int kbe  = kb;
-      unsigned int ktw  = kw;
-      unsigned int kbs  = neighborZ[ks];
-      unsigned int ktn  = k;
-      unsigned int kbn  = kb;
-      unsigned int kts  = ks;
-      unsigned int ktse = ks;
-      unsigned int kbnw = kbw;
-      unsigned int ktnw = kw;
-      unsigned int kbse = kbs;
-      unsigned int ktsw = ksw;
-      unsigned int kbne = kb;
-      unsigned int ktne = k;
-      unsigned int kbsw = neighborZ[ksw];
-      //////////////////////////////////////////////////////////////////////////
-      pressD[k] = c0o1;
-	  rhoD[k]   = c0o1;
-	  vxD[k]    = c0o1;
-	  vyD[k]    = c0o1;
-	  vzD[k]    = c0o1;
-
-      if(geoD[k] == GEO_FLUID)
-      {
-         rhoD[k]    =   (D.f[dirE   ])[ke  ]+ (D.f[dirW   ])[kw  ]+ 
-                        (D.f[dirN   ])[kn  ]+ (D.f[dirS   ])[ks  ]+
-                        (D.f[dirT   ])[kt  ]+ (D.f[dirB   ])[kb  ]+
-                        (D.f[dirNE  ])[kne ]+ (D.f[dirSW  ])[ksw ]+
-                        (D.f[dirSE  ])[kse ]+ (D.f[dirNW  ])[knw ]+
-                        (D.f[dirTE  ])[kte ]+ (D.f[dirBW  ])[kbw ]+
-                        (D.f[dirBE  ])[kbe ]+ (D.f[dirTW  ])[ktw ]+
-                        (D.f[dirTN  ])[ktn ]+ (D.f[dirBS  ])[kbs ]+
-                        (D.f[dirBN  ])[kbn ]+ (D.f[dirTS  ])[kts ]+
-                        (D.f[dirZERO])[kzero]+ 
-                        (D.f[dirTNE ])[ktne]+ (D.f[dirTSW ])[ktsw]+ 
-                        (D.f[dirTSE ])[ktse]+ (D.f[dirTNW ])[ktnw]+ 
-                        (D.f[dirBNE ])[kbne]+ (D.f[dirBSW ])[kbsw]+ 
-                        (D.f[dirBSE ])[kbse]+ (D.f[dirBNW ])[kbnw];
-
-         vxD[k]     =  ((D.f[dirE   ])[ke  ]- (D.f[dirW   ])[kw  ]+ 
-                        (D.f[dirNE  ])[kne ]- (D.f[dirSW  ])[ksw ]+
-                        (D.f[dirSE  ])[kse ]- (D.f[dirNW  ])[knw ]+
-                        (D.f[dirTE  ])[kte ]- (D.f[dirBW  ])[kbw ]+
-                        (D.f[dirBE  ])[kbe ]- (D.f[dirTW  ])[ktw ]+
-                        (D.f[dirTNE ])[ktne]- (D.f[dirTSW ])[ktsw]+ 
-                        (D.f[dirTSE ])[ktse]- (D.f[dirTNW ])[ktnw]+ 
-                        (D.f[dirBNE ])[kbne]- (D.f[dirBSW ])[kbsw]+ 
-                        (D.f[dirBSE ])[kbse]- (D.f[dirBNW ])[kbnw])/(c1o1+rhoD[k]);
-
-         vyD[k]     =  ((D.f[dirN   ])[kn  ]- (D.f[dirS   ])[ks  ]+
-                        (D.f[dirNE  ])[kne ]- (D.f[dirSW  ])[ksw ]-
-                        (D.f[dirSE  ])[kse ]+ (D.f[dirNW  ])[knw ]+
-                        (D.f[dirTN  ])[ktn ]- (D.f[dirBS  ])[kbs ]+
-                        (D.f[dirBN  ])[kbn ]- (D.f[dirTS  ])[kts ]+
-                        (D.f[dirTNE ])[ktne]- (D.f[dirTSW ])[ktsw]- 
-                        (D.f[dirTSE ])[ktse]+ (D.f[dirTNW ])[ktnw]+ 
-                        (D.f[dirBNE ])[kbne]- (D.f[dirBSW ])[kbsw]- 
-                        (D.f[dirBSE ])[kbse]+ (D.f[dirBNW ])[kbnw])/(c1o1+rhoD[k]);
-
-         vzD[k]     =  ((D.f[dirT   ])[kt  ]- (D.f[dirB   ])[kb  ]+
-                        (D.f[dirTE  ])[kte ]- (D.f[dirBW  ])[kbw ]-
-                        (D.f[dirBE  ])[kbe ]+ (D.f[dirTW  ])[ktw ]+
-                        (D.f[dirTN  ])[ktn ]- (D.f[dirBS  ])[kbs ]-
-                        (D.f[dirBN  ])[kbn ]+ (D.f[dirTS  ])[kts ]+
-                        (D.f[dirTNE ])[ktne]+ (D.f[dirTSW ])[ktsw]+ 
-                        (D.f[dirTSE ])[ktse]+ (D.f[dirTNW ])[ktnw]- 
-                        (D.f[dirBNE ])[kbne]- (D.f[dirBSW ])[kbsw]- 
-                        (D.f[dirBSE ])[kbse]- (D.f[dirBNW ])[kbnw])/(c1o1+rhoD[k]);
-
-         pressD[k]  =  ((D.f[dirE   ])[ke  ]+ (D.f[dirW   ])[kw  ]+ 
-                        (D.f[dirN   ])[kn  ]+ (D.f[dirS   ])[ks  ]+
-                        (D.f[dirT   ])[kt  ]+ (D.f[dirB   ])[kb  ]+
-                        2.f*(
-                        (D.f[dirNE  ])[kne ]+ (D.f[dirSW  ])[ksw ]+
-                        (D.f[dirSE  ])[kse ]+ (D.f[dirNW  ])[knw ]+
-                        (D.f[dirTE  ])[kte ]+ (D.f[dirBW  ])[kbw ]+
-                        (D.f[dirBE  ])[kbe ]+ (D.f[dirTW  ])[ktw ]+
-                        (D.f[dirTN  ])[ktn ]+ (D.f[dirBS  ])[kbs ]+
-                        (D.f[dirBN  ])[kbn ]+ (D.f[dirTS  ])[kts ])+
-                        3.f*(
-                        (D.f[dirTNE ])[ktne]+ (D.f[dirTSW ])[ktsw]+ 
-                        (D.f[dirTSE ])[ktse]+ (D.f[dirTNW ])[ktnw]+ 
-                        (D.f[dirBNE ])[kbne]+ (D.f[dirBSW ])[kbsw]+ 
-                        (D.f[dirBSE ])[kbse]+ (D.f[dirBNW ])[kbnw])-
-                        rhoD[k]-(vxD[k] * vxD[k] + vyD[k] * vyD[k] + vzD[k] * vzD[k]) * (c1o1+rhoD[k])) * c1o2+rhoD[k]; // times zero for incompressible case   
-         //achtung op hart gesetzt Annahme op = 1 ;                                                    ^^^^(1.0/op-0.5)=0.5
-
-      }
-   }
-
+    rhoD[k]   = vf::lbm::getDensity(distribution.f);
+    vxD[k]    = vf::lbm::getCompressibleVelocityX1(distribution.f, rhoD[k]);
+    vyD[k]    = vf::lbm::getCompressibleVelocityX2(distribution.f, rhoD[k]);
+    vzD[k]    = vf::lbm::getCompressibleVelocityX3(distribution.f, rhoD[k]);
+    pressD[k] = vf::lbm::getPressure(distribution.f, rhoD[k], vxD[k], vyD[k], vzD[k]); 
 }
 
 ////////////////////////////////////////////////////////////////////////////////
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/RunLBMKernel.cuh b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/RunLBMKernel.cuh
index b4097851b251b7447f6ce06856d0b9187999a20b..ecfdb8e38bc5da02a473b6f2e67340963ed8b38a 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/RunLBMKernel.cuh
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Kernels/BasicKernels/FluidFlow/Compressible/RunLBMKernel.cuh
@@ -32,9 +32,11 @@ template<typename KernelFunctor>
 __global__ void runKernel(KernelFunctor kernel, GPUKernelParameter kernelParameter)
 {
     const uint k = getNodeIndex();
-    const uint nodeType = kernelParameter.typeOfGridNode[k];
 
-    if (!isValidFluidNode(k, kernelParameter.size_Mat, nodeType))
+    if(k >= kernelParameter.size_Mat)
+        return;
+
+    if (!isValidFluidNode(kernelParameter.typeOfGridNode[k]))
         return;
 
     DistributionWrapper distributionWrapper {
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cu b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cu
index bbb01d95410612d36d55f1e0113175a8741b9ade..2b9203bf96be89fb328cbbe221994dffed481e14 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cu
@@ -158,10 +158,9 @@ __device__ unsigned int getNodeIndex()
     return nx * (ny * z + y) + x;
 }
 
-__device__ bool isValidFluidNode(uint k, int size_Mat, uint nodeType)
+__device__ bool isValidFluidNode(uint nodeType)
 {
-    return (k < size_Mat) &&
-           (nodeType == GEO_FLUID || nodeType == GEO_PM_0 || nodeType == GEO_PM_1 || nodeType == GEO_PM_2);
+    return (nodeType == GEO_FLUID || nodeType == GEO_PM_0 || nodeType == GEO_PM_1 || nodeType == GEO_PM_2);
 }
 
 
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cuh b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cuh
index 935030701924233d959fb69b74a7c3087feb0834..6b38cac75c99680c71420533455dd060195b6c87 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cuh
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cuh
@@ -90,7 +90,7 @@ struct DistributionWrapper
 
 __device__ unsigned int getNodeIndex();
 
-__device__ bool isValidFluidNode(uint k, int size_Mat, uint nodeType);
+__device__ bool isValidFluidNode(uint nodeType);
 
 }
 }