diff --git a/src/gpu/core/BoundaryConditions/BoundaryConditionFactory.cpp b/src/gpu/core/BoundaryConditions/BoundaryConditionFactory.cpp
index b1c398638cff1ec1b6d52f59f8e773183e270331..6e80f71fe45abeddce4a11b9dedd2b476b595286 100644
--- a/src/gpu/core/BoundaryConditions/BoundaryConditionFactory.cpp
+++ b/src/gpu/core/BoundaryConditions/BoundaryConditionFactory.cpp
@@ -1,9 +1,12 @@
 #include "BoundaryConditionFactory.h"
-#include "GPU/GPU_Interface.h"
 #include "Parameter/Parameter.h"
 #include "grid/BoundaryConditions/BoundaryCondition.h"
 #include <variant>
 
+#include "BoundaryConditions/Outflow/Outflow.h"
+#include "GPU/GPU_Interface.h"
+
+
 void BoundaryConditionFactory::setVelocityBoundaryCondition(VelocityBC boundaryConditionType)
 {
     this->velocityBoundaryCondition = boundaryConditionType;
@@ -135,10 +138,10 @@ boundaryCondition BoundaryConditionFactory::getPressureBoundaryConditionPre() co
             return QPressDevNEQ27;
             break;
         case PressureBC::OutflowNonReflective:
-            return QPressNoRhoDev27;
+            return OutflowNonReflecting;
             break;
         case PressureBC::OutflowNonReflectivePressureCorrection:
-            return QPressZeroRhoOutflowDev27;
+            return OutflowNonReflectingPressureCorrection;
         default:
             return nullptr;
     }
diff --git a/src/gpu/core/BoundaryConditions/BoundaryConditionFactoryTest.cpp b/src/gpu/core/BoundaryConditions/BoundaryConditionFactoryTest.cpp
index a01b6af7abfebcdf488a66fc0c0ba1c9d3c81987..024a5fad27453a3c61243a3e3067acc2342e3669 100644
--- a/src/gpu/core/BoundaryConditions/BoundaryConditionFactoryTest.cpp
+++ b/src/gpu/core/BoundaryConditions/BoundaryConditionFactoryTest.cpp
@@ -2,9 +2,11 @@
 #include <typeindex>
 
 #include "BoundaryConditionFactory.h"
-#include "GPU/GPU_Interface.h"
 #include "gpu/GridGenerator/grid/BoundaryConditions/BoundaryCondition.h"
 
+#include "BoundaryConditions/Outflow/Outflow.h"
+#include "GPU/GPU_Interface.h"
+
 using bcFunction = void (*)(LBMSimulationParameter *, QforBoundaryConditions *);
 using bcFunctionParamter = void (*)(Parameter *, QforBoundaryConditions *, const int level);
 
@@ -173,8 +175,8 @@ TEST(BoundaryConditionFactoryTest, pressureBC)
         << "The returned boundary condition is not the expected function QPressDevNEQ27.";
 
     bcFactory.setPressureBoundaryCondition(BoundaryConditionFactory::PressureBC::OutflowNonReflective);
-    EXPECT_TRUE( *(getPressureBcTarget(bcFactory)) == QPressNoRhoDev27)
-        << "The returned boundary condition is not the expected function QPressNoRhoDev27.";
+    EXPECT_TRUE(*(getPressureBcTarget(bcFactory)) == OutflowNonReflecting)
+        << "The returned boundary condition is not the expected function OutflowNonReflecting_Device.";
 }
 
 bcFunction getGeometryBcTarget(BoundaryConditionFactory &bcFactory)
diff --git a/src/gpu/core/BoundaryConditions/BoundaryConditionKernelManager.cpp b/src/gpu/core/BoundaryConditions/BoundaryConditionKernelManager.cpp
index 60cd1f1138f7aa9362c07332f7307cd13d5ce3c1..4aa4983c08dffa2b2d29138cf9bed4dd4ac589db 100644
--- a/src/gpu/core/BoundaryConditions/BoundaryConditionKernelManager.cpp
+++ b/src/gpu/core/BoundaryConditions/BoundaryConditionKernelManager.cpp
@@ -41,9 +41,11 @@
 #include "GridGenerator/TransientBCSetter/TransientBCSetter.h"
 #include "Calculation/Cp.h"
 #include "Calculation/DragLift.h"
-#include "GPU/GPU_Interface.h"
 #include "Parameter/Parameter.h"
 
+#include "BoundaryConditions/Outflow/Outflow.h"
+#include "GPU/GPU_Interface.h"
+
 BoundaryConditionKernelManager::BoundaryConditionKernelManager(SPtr<Parameter> parameter, BoundaryConditionFactory *bcFactory) : para(parameter)
 {
     this->velocityBoundaryConditionPost = bcFactory->getVelocityBoundaryConditionPost();
@@ -321,7 +323,7 @@ void BoundaryConditionKernelManager::runGeoBCKernelPost(const int level) const
 void BoundaryConditionKernelManager::runOutflowBCKernelPre(const int level) const{
     if (para->getParD(level)->outflowBC.numberOfBCnodes > 0)
     {
-        QPressNoRhoDev27(para->getParD(level).get(), &(para->getParD(level)->outflowBC));
+        OutflowNonReflecting(para->getParD(level).get(), &(para->getParD(level)->outflowBC));
 
         // TODO: https://git.rz.tu-bs.de/irmb/VirtualFluids_dev/-/issues/29
         // if (  myid == numprocs - 1)
diff --git a/src/gpu/core/BoundaryConditions/Outflow/Outflow.cu b/src/gpu/core/BoundaryConditions/Outflow/Outflow.cu
new file mode 100644
index 0000000000000000000000000000000000000000..4a492e93ec07bdf89d79588c40cedc5cd47272e4
--- /dev/null
+++ b/src/gpu/core/BoundaryConditions/Outflow/Outflow.cu
@@ -0,0 +1,59 @@
+//  _    ___      __              __________      _     __        ______________   __
+// | |  / (_)____/ /___  ______ _/ / ____/ /_  __(_)___/ /____   /  ___/ __  / /  / /
+// | | / / / ___/ __/ / / / __ `/ / /_  / / / / / / __  / ___/  / /___/ /_/ / /  / /
+// | |/ / / /  / /_/ /_/ / /_/ / / __/ / / /_/ / / /_/ (__  )  / /_) / ____/ /__/ /
+// |___/_/_/   \__/\__,_/\__,_/_/_/   /_/\__,_/_/\__,_/____/   \____/_/    \_____/
+//
+//////////////////////////////////////////////////////////////////////////
+#include <cuda_runtime.h>
+#include <helper_functions.h>
+#include <helper_cuda.h>
+
+#include "LBM/LB.h"
+#include <cuda_helper/CudaGrid.h>
+
+#include "BoundaryConditions/Outflow/Outflow_Device.cuh"
+#include "Parameter/Parameter.h"
+
+void OutflowNonReflecting(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition)
+{
+    dim3 grid = vf::cuda::getCudaGrid( parameterDevice->numberofthreads,  boundaryCondition->numberOfBCnodes);
+    dim3 threads(parameterDevice->numberofthreads, 1, 1 );
+
+    OutflowNonReflecting_Device<<< grid, threads >>> (
+        boundaryCondition->RhoBC,
+        parameterDevice->distributions.f[0],
+        boundaryCondition->k,
+        boundaryCondition->kN,
+        boundaryCondition->numberOfBCnodes,
+        parameterDevice->omega,
+        parameterDevice->neighborX,
+        parameterDevice->neighborY,
+        parameterDevice->neighborZ,
+        parameterDevice->numberOfNodes,
+        parameterDevice->isEvenTimestep,
+        vf::lbm::dir::dP00);
+    getLastCudaError("OutflowNonReflecting_Device execution failed");
+}
+
+void OutflowNonReflectingPressureCorrection(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition)
+{
+    dim3 grid = vf::cuda::getCudaGrid( parameterDevice->numberofthreads,  boundaryCondition->numberOfBCnodes);
+    dim3 threads(parameterDevice->numberofthreads, 1, 1 );
+
+    OutflowNonReflectingPressureCorrection_Device<<< grid, threads >>> (
+        boundaryCondition->RhoBC,
+        parameterDevice->distributions.f[0],
+        boundaryCondition->k,
+        boundaryCondition->kN,
+        boundaryCondition->numberOfBCnodes,
+        parameterDevice->omega,
+        parameterDevice->neighborX,
+        parameterDevice->neighborY,
+        parameterDevice->neighborZ,
+        parameterDevice->numberOfNodes,
+        parameterDevice->isEvenTimestep,
+        vf::lbm::dir::dP00,
+        parameterDevice->outflowPressureCorrectionFactor);
+    getLastCudaError("OutflowNonReflectingPressureCorrection_Device execution failed");
+}
diff --git a/src/gpu/core/BoundaryConditions/Outflow/Outflow.h b/src/gpu/core/BoundaryConditions/Outflow/Outflow.h
new file mode 100644
index 0000000000000000000000000000000000000000..5a885d1dd4b9959cad75e250ae66f46a452ce7e6
--- /dev/null
+++ b/src/gpu/core/BoundaryConditions/Outflow/Outflow.h
@@ -0,0 +1,19 @@
+//  _    ___      __              __________      _     __        ______________   __
+// | |  / (_)____/ /___  ______ _/ / ____/ /_  __(_)___/ /____   /  ___/ __  / /  / /
+// | | / / / ___/ __/ / / / __ `/ / /_  / / / / / / __  / ___/  / /___/ /_/ / /  / /
+// | |/ / / /  / /_/ /_/ / /_/ / / __/ / / /_/ / / /_/ (__  )  / /_) / ____/ /__/ / 
+// |___/_/_/   \__/\__,_/\__,_/_/_/   /_/\__,_/_/\__,_/____/   \____/_/    \_____/
+//
+//////////////////////////////////////////////////////////////////////////
+#ifndef Outflow_H
+#define Outflow_H
+
+#include "LBM/LB.h"
+
+struct LBMSimulationParameter;
+
+void OutflowNonReflecting(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition);
+
+void OutflowNonReflectingPressureCorrection(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition);
+
+#endif
diff --git a/src/gpu/core/BoundaryConditions/Outflow/OutflowNonReflecting.cu b/src/gpu/core/BoundaryConditions/Outflow/OutflowNonReflecting.cu
new file mode 100644
index 0000000000000000000000000000000000000000..11a891dbc29699748d9b132f8481dc97a1dd3c73
--- /dev/null
+++ b/src/gpu/core/BoundaryConditions/Outflow/OutflowNonReflecting.cu
@@ -0,0 +1,274 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+//      \    \  |    |   ________________________________________________________________
+//       \    \ |    |  |  ______________________________________________________________|
+//        \    \|    |  |  |         __          __     __     __     ______      _______
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \author Martin Schoenherr, Anna Wellmann
+//======================================================================================
+#include "LBM/LB.h"
+#include "lbm/constants/D3Q27.h"
+#include "basics/constants/NumericConstants.h"
+#include "LBM/GPUHelperFunctions/KernelUtilities.h"
+
+using namespace vf::basics::constant;
+using namespace vf::lbm::dir;
+using namespace vf::gpu;
+
+__host__ __device__ real computeOutflowDistribution(const real* const &f, const real* const &f1, const int dir, const real cs)
+{
+   return f1[dir] * cs + (c1o1 - cs) * f[dir];
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+__global__ void OutflowNonReflecting_Device(
+    real* rhoBC,
+    real* distributions,
+    int* k_Q,
+    int* k_N,
+    int numberOfBCnodes,
+    real om1,
+    unsigned int* neighborX,
+    unsigned int* neighborY,
+    unsigned int* neighborZ,
+    unsigned long long numberOfLBnodes,
+    bool isEvenTimestep,
+    int direction)
+{
+   ////////////////////////////////////////////////////////////////////////////////
+   //! - Get the node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
+   //!
+   const unsigned nodeIndex = getNodeIndex();
+
+   //////////////////////////////////////////////////////////////////////////
+
+   if(nodeIndex >= numberOfBCnodes) return;
+
+   ////////////////////////////////////////////////////////////////////////////////
+   //index
+   unsigned int KQK  = k_Q[nodeIndex];
+   // unsigned int kzero= KQK;
+   unsigned int ke   = KQK;
+   unsigned int kw   = neighborX[KQK];
+   unsigned int kn   = KQK;
+   unsigned int ks   = neighborY[KQK];
+   unsigned int kt   = KQK;
+   unsigned int kb   = neighborZ[KQK];
+   unsigned int ksw  = neighborY[kw];
+   unsigned int kne  = KQK;
+   unsigned int kse  = ks;
+   unsigned int knw  = kw;
+   unsigned int kbw  = neighborZ[kw];
+   unsigned int kte  = KQK;
+   unsigned int kbe  = kb;
+   unsigned int ktw  = kw;
+   unsigned int kbs  = neighborZ[ks];
+   unsigned int ktn  = KQK;
+   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 = KQK;
+   unsigned int kbsw = neighborZ[ksw];
+   ////////////////////////////////////////////////////////////////////////////////
+   //index1
+   unsigned int K1QK  = k_N[nodeIndex];
+   //unsigned int k1zero= K1QK;
+   unsigned int k1e   = K1QK;
+   unsigned int k1w   = neighborX[K1QK];
+   unsigned int k1n   = K1QK;
+   unsigned int k1s   = neighborY[K1QK];
+   unsigned int k1t   = K1QK;
+   unsigned int k1b   = neighborZ[K1QK];
+   unsigned int k1sw  = neighborY[k1w];
+   unsigned int k1ne  = K1QK;
+   unsigned int k1se  = k1s;
+   unsigned int k1nw  = k1w;
+   unsigned int k1bw  = neighborZ[k1w];
+   unsigned int k1te  = K1QK;
+   unsigned int k1be  = k1b;
+   unsigned int k1tw  = k1w;
+   unsigned int k1bs  = neighborZ[k1s];
+   unsigned int k1tn  = K1QK;
+   unsigned int k1bn  = k1b;
+   unsigned int k1ts  = k1s;
+   unsigned int k1tse = k1s;
+   unsigned int k1bnw = k1bw;
+   unsigned int k1tnw = k1w;
+   unsigned int k1bse = k1bs;
+   unsigned int k1tsw = k1sw;
+   unsigned int k1bne = k1b;
+   unsigned int k1tne = K1QK;
+   unsigned int k1bsw = neighborZ[k1sw];
+   ////////////////////////////////////////////////////////////////////////////////
+   Distributions27 dist;
+   getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
+   real f[27], f1[27];
+   //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+   f1[dP00] = (dist.f[dP00])[k1e   ];
+   f1[dM00] = (dist.f[dM00])[k1w   ];
+   f1[d0P0] = (dist.f[d0P0])[k1n   ];
+   f1[d0M0] = (dist.f[d0M0])[k1s   ];
+   f1[d00P] = (dist.f[d00P])[k1t   ];
+   f1[d00M] = (dist.f[d00M])[k1b   ];
+   f1[dPP0] = (dist.f[dPP0])[k1ne  ];
+   f1[dMM0] = (dist.f[dMM0])[k1sw  ];
+   f1[dPM0] = (dist.f[dPM0])[k1se  ];
+   f1[dMP0] = (dist.f[dMP0])[k1nw  ];
+   f1[dP0P] = (dist.f[dP0P])[k1te  ];
+   f1[dM0M] = (dist.f[dM0M])[k1bw  ];
+   f1[dP0M] = (dist.f[dP0M])[k1be  ];
+   f1[dM0P] = (dist.f[dM0P])[k1tw  ];
+   f1[d0PP] = (dist.f[d0PP])[k1tn  ];
+   f1[d0MM] = (dist.f[d0MM])[k1bs  ];
+   f1[d0PM] = (dist.f[d0PM])[k1bn  ];
+   f1[d0MP] = (dist.f[d0MP])[k1ts  ];
+   // f1[d000] = (dist.f[d000])[k1zero];
+   f1[dPPP] = (dist.f[dPPP])[k1tne ];
+   f1[dMMP] = (dist.f[dMMP])[k1tsw ];
+   f1[dPMP] = (dist.f[dPMP])[k1tse ];
+   f1[dMPP] = (dist.f[dMPP])[k1tnw ];
+   f1[dPPM] = (dist.f[dPPM])[k1bne ];
+   f1[dMMM] = (dist.f[dMMM])[k1bsw ];
+   f1[dPMM] = (dist.f[dPMM])[k1bse ];
+   f1[dMPM] = (dist.f[dMPM])[k1bnw ];
+   //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+   f[dP00] = (dist.f[dP00])[ke   ];
+   f[dM00] = (dist.f[dM00])[kw   ];
+   f[d0P0] = (dist.f[d0P0])[kn   ];
+   f[d0M0] = (dist.f[d0M0])[ks   ];
+   f[d00P] = (dist.f[d00P])[kt   ];
+   f[d00M] = (dist.f[d00M])[kb   ];
+   f[dPP0] = (dist.f[dPP0])[kne  ];
+   f[dMM0] = (dist.f[dMM0])[ksw  ];
+   f[dPM0] = (dist.f[dPM0])[kse  ];
+   f[dMP0] = (dist.f[dMP0])[knw  ];
+   f[dP0P] = (dist.f[dP0P])[kte  ];
+   f[dM0M] = (dist.f[dM0M])[kbw  ];
+   f[dP0M] = (dist.f[dP0M])[kbe  ];
+   f[dM0P] = (dist.f[dM0P])[ktw  ];
+   f[d0PP] = (dist.f[d0PP])[ktn  ];
+   f[d0MM] = (dist.f[d0MM])[kbs  ];
+   f[d0PM] = (dist.f[d0PM])[kbn  ];
+   f[d0MP] = (dist.f[d0MP])[kts  ];
+   // f[d000] = (dist.f[d000])[kzero];
+   f[dPPP] = (dist.f[dPPP])[ktne ];
+   f[dMMP] = (dist.f[dMMP])[ktsw ];
+   f[dPMP] = (dist.f[dPMP])[ktse ];
+   f[dMPP] = (dist.f[dMPP])[ktnw ];
+   f[dPPM] = (dist.f[dPPM])[kbne ];
+   f[dMMM] = (dist.f[dMMM])[kbsw ];
+   f[dPMM] = (dist.f[dPMM])[kbse ];
+   f[dMPM] = (dist.f[dMPM])[kbnw ];
+   //////////////////////////////////////////////////////////////////////////
+
+
+   real cs = c1o1 / sqrtf(c3o1);
+
+   //////////////////////////////////////////////////////////////////////////
+   getPointersToDistributions(dist, distributions, numberOfLBnodes, !isEvenTimestep);
+   switch(direction)
+   {
+      case dM00:
+         (dist.f[dP00])[ke   ] = computeOutflowDistribution(f, f1, dP00, cs);
+         (dist.f[dPM0])[kse  ] = computeOutflowDistribution(f, f1, dPM0, cs);
+         (dist.f[dPP0])[kne  ] = computeOutflowDistribution(f, f1, dPP0, cs);
+         (dist.f[dP0M])[kbe  ] = computeOutflowDistribution(f, f1, dP0M, cs);
+         (dist.f[dP0P])[kte  ] = computeOutflowDistribution(f, f1, dP0P, cs);
+         (dist.f[dPMP])[ktse ] = computeOutflowDistribution(f, f1, dPMP, cs);
+         (dist.f[dPPP])[ktne ] = computeOutflowDistribution(f, f1, dPPP, cs);
+         (dist.f[dPMM])[kbse ] = computeOutflowDistribution(f, f1, dPMM, cs);
+         (dist.f[dPPM])[kbne ] = computeOutflowDistribution(f, f1, dPPM, cs);
+         break;
+
+      case dP00:
+         (dist.f[dM00])[kw   ] = computeOutflowDistribution(f, f1, dM00, cs);
+         (dist.f[dMM0])[ksw  ] = computeOutflowDistribution(f, f1, dMM0, cs);
+         (dist.f[dMP0])[knw  ] = computeOutflowDistribution(f, f1, dMP0, cs);
+         (dist.f[dM0M])[kbw  ] = computeOutflowDistribution(f, f1, dM0M, cs);
+         (dist.f[dM0P])[ktw  ] = computeOutflowDistribution(f, f1, dM0P, cs);
+         (dist.f[dMMP])[ktsw ] = computeOutflowDistribution(f, f1, dMMP, cs);
+         (dist.f[dMPP])[ktnw ] = computeOutflowDistribution(f, f1, dMPP, cs);
+         (dist.f[dMMM])[kbsw ] = computeOutflowDistribution(f, f1, dMMM, cs);
+         (dist.f[dMPM])[kbnw ] = computeOutflowDistribution(f, f1, dMPM, cs);
+         break;
+
+      case d0M0:
+         (dist.f[d0P0])[kn   ] = computeOutflowDistribution(f, f1, d0P0, cs);
+         (dist.f[dPP0])[kne  ] = computeOutflowDistribution(f, f1, dPP0, cs);
+         (dist.f[dMP0])[knw  ] = computeOutflowDistribution(f, f1, dMP0, cs);
+         (dist.f[d0PP])[ktn  ] = computeOutflowDistribution(f, f1, d0PP, cs);
+         (dist.f[d0PM])[kbn  ] = computeOutflowDistribution(f, f1, d0PM, cs);
+         (dist.f[dPPP])[ktne ] = computeOutflowDistribution(f, f1, dPPP, cs);
+         (dist.f[dMPP])[ktnw ] = computeOutflowDistribution(f, f1, dMPP, cs);
+         (dist.f[dPPM])[kbne ] = computeOutflowDistribution(f, f1, dPPM, cs);
+         (dist.f[dMPM])[kbnw ] = computeOutflowDistribution(f, f1, dMPM, cs);
+         break;
+
+      case d0P0:
+         (dist.f[d0M0])[ks   ] = computeOutflowDistribution(f, f1, d0M0, cs);
+         (dist.f[dPM0])[kse  ] = computeOutflowDistribution(f, f1, dPM0, cs);
+         (dist.f[dMM0])[ksw  ] = computeOutflowDistribution(f, f1, dMM0, cs);
+         (dist.f[d0MP])[kts  ] = computeOutflowDistribution(f, f1, d0MP, cs);
+         (dist.f[d0MM])[kbs  ] = computeOutflowDistribution(f, f1, d0MM, cs);
+         (dist.f[dPMP])[ktse ] = computeOutflowDistribution(f, f1, dPMP, cs);
+         (dist.f[dMMP])[ktsw ] = computeOutflowDistribution(f, f1, dMMP, cs);
+         (dist.f[dPMM])[kbse ] = computeOutflowDistribution(f, f1, dPMM, cs);
+         (dist.f[dMMM])[kbsw ] = computeOutflowDistribution(f, f1, dMMM, cs);
+         break;
+
+      case d00M:
+         (dist.f[d00P])[kt   ] = computeOutflowDistribution(f, f1, d00P, cs);
+         (dist.f[dP0P])[kte  ] = computeOutflowDistribution(f, f1, dP0P, cs);
+         (dist.f[dM0P])[ktw  ] = computeOutflowDistribution(f, f1, dM0P, cs);
+         (dist.f[d0PP])[ktn  ] = computeOutflowDistribution(f, f1, d0PP, cs);
+         (dist.f[d0MP])[kts  ] = computeOutflowDistribution(f, f1, d0MP, cs);
+         (dist.f[dPPP])[ktne ] = computeOutflowDistribution(f, f1, dPPP, cs);
+         (dist.f[dMPP])[ktnw ] = computeOutflowDistribution(f, f1, dMPP, cs);
+         (dist.f[dPMP])[ktse ] = computeOutflowDistribution(f, f1, dPMP, cs);
+         (dist.f[dMMP])[ktsw ] = computeOutflowDistribution(f, f1, dMMP, cs);
+         break;
+
+      case d00P:
+         (dist.f[d00M])[kb   ] = computeOutflowDistribution(f, f1, d00M, cs);
+         (dist.f[dP0M])[kbe  ] = computeOutflowDistribution(f, f1, dP0M, cs);
+         (dist.f[dM0M])[kbw  ] = computeOutflowDistribution(f, f1, dM0M, cs);
+         (dist.f[d0PM])[kbn  ] = computeOutflowDistribution(f, f1, d0PM, cs);
+         (dist.f[d0MM])[kbs  ] = computeOutflowDistribution(f, f1, d0MM, cs);
+         (dist.f[dPPM])[kbne ] = computeOutflowDistribution(f, f1, dPPM, cs);
+         (dist.f[dMPM])[kbnw ] = computeOutflowDistribution(f, f1, dMPM, cs);
+         (dist.f[dPMM])[kbse ] = computeOutflowDistribution(f, f1, dPMM, cs);
+         (dist.f[dMMM])[kbsw ] = computeOutflowDistribution(f, f1, dMMM, cs);
+         break;
+      default:
+         break;
+   }
+}
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
diff --git a/src/gpu/core/BoundaryConditions/Outflow/OutflowNonReflectingPressureCorrection.cu b/src/gpu/core/BoundaryConditions/Outflow/OutflowNonReflectingPressureCorrection.cu
new file mode 100644
index 0000000000000000000000000000000000000000..e70d4721d954da4ca95b140e931b789a2f3cf7d4
--- /dev/null
+++ b/src/gpu/core/BoundaryConditions/Outflow/OutflowNonReflectingPressureCorrection.cu
@@ -0,0 +1,239 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+//      \    \  |    |   ________________________________________________________________
+//       \    \ |    |  |  ______________________________________________________________|
+//        \    \|    |  |  |         __          __     __     __     ______      _______
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \author Martin Schoenherr, Anna Wellmann
+//======================================================================================
+#include "LBM/LB.h"
+#include "lbm/constants/D3Q27.h"
+#include "basics/constants/NumericConstants.h"
+#include "lbm/MacroscopicQuantities.h"
+#include "LBM/GPUHelperFunctions/KernelUtilities.h"
+
+using namespace vf::basics::constant;
+using namespace vf::lbm::dir;
+using namespace vf::gpu;
+
+
+__host__ __device__ real computeOutflowDistribution(const real* const &f, const real* const &f1, const int dir, const real rhoCorrection, const real cs, const real weight)
+{
+   return f1[dir  ] * cs + (c1o1 - cs) * f[dir  ] - weight *rhoCorrection;
+}
+
+__global__ void OutflowNonReflectingPressureCorrection_Device(
+    real* rhoBC,
+    real* distributions,
+    int* k_Q,
+    int* k_N,
+    int numberOfBCnodes,
+    real om1,
+    unsigned int* neighborX,
+    unsigned int* neighborY,
+    unsigned int* neighborZ,
+    unsigned long long numberOfLBnodes,
+    bool isEvenTimestep,
+    int direction,
+    real densityCorrectionFactor)
+{
+   ////////////////////////////////////////////////////////////////////////////////
+   //! - Get the node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
+   //!
+   const unsigned nodeIndex = getNodeIndex();
+
+   //////////////////////////////////////////////////////////////////////////
+
+   if( nodeIndex >= numberOfBCnodes ) return;
+
+   ////////////////////////////////////////////////////////////////////////////////
+   //index
+
+   uint k_000 = k_Q[nodeIndex];
+   uint k_M00 = neighborX[k_000];
+   uint k_0M0 = neighborY[k_000];
+   uint k_00M = neighborZ[k_000];
+   uint k_MM0 = neighborY[k_M00];
+   uint k_M0M = neighborZ[k_M00];
+   uint k_0MM = neighborZ[k_0M0];
+   uint k_MMM = neighborZ[k_MM0];
+
+   ////////////////////////////////////////////////////////////////////////////////
+   //index of neighbor
+   uint kN_000 = k_N[nodeIndex];
+   uint kN_M00 = neighborX[k_000];
+   uint kN_0M0 = neighborY[k_000];
+   uint kN_00M = neighborZ[k_000];
+   uint kN_MM0 = neighborY[k_M00];
+   uint kN_M0M = neighborZ[k_M00];
+   uint kN_0MM = neighborZ[k_0M0];
+   uint kN_MMM = neighborZ[k_MM0];
+   ////////////////////////////////////////////////////////////////////////////////
+   Distributions27 dist;
+   getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
+   real f[27], fN[27];
+   //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+   f[d000] = (dist.f[d000])[k_000];
+   f[dP00] = (dist.f[dP00])[k_000];
+   f[dM00] = (dist.f[dM00])[k_M00];
+   f[d0P0] = (dist.f[d0P0])[k_000];
+   f[d0M0] = (dist.f[d0M0])[k_0M0];
+   f[d00P] = (dist.f[d00P])[k_000];
+   f[d00M] = (dist.f[d00M])[k_00M];
+   f[dPP0] = (dist.f[dPP0])[k_000];
+   f[dMM0] = (dist.f[dMM0])[k_MM0];
+   f[dPM0] = (dist.f[dPM0])[k_0M0];
+   f[dMP0] = (dist.f[dMP0])[k_M00];
+   f[dP0P] = (dist.f[dP0P])[k_000];
+   f[dM0M] = (dist.f[dM0M])[k_M0M];
+   f[dP0M] = (dist.f[dP0M])[k_00M];
+   f[dM0P] = (dist.f[dM0P])[k_M00];
+   f[d0PP] = (dist.f[d0PP])[k_000];
+   f[d0MM] = (dist.f[d0MM])[k_0MM];
+   f[d0PM] = (dist.f[d0PM])[k_00M];
+   f[d0MP] = (dist.f[d0MP])[k_0M0];
+   f[dPPP] = (dist.f[dPPP])[k_000];
+   f[dMPP] = (dist.f[dMPP])[k_M00];
+   f[dPMP] = (dist.f[dPMP])[k_0M0];
+   f[dMMP] = (dist.f[dMMP])[k_MM0];
+   f[dPPM] = (dist.f[dPPM])[k_00M];
+   f[dMPM] = (dist.f[dMPM])[k_M0M];
+   f[dPMM] = (dist.f[dPMM])[k_0MM];
+   f[dMMM] = (dist.f[dMMM])[k_MMM];
+   //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+   fN[d000] = (dist.f[d000])[kN_000];
+   fN[dP00] = (dist.f[dP00])[kN_000];
+   fN[dM00] = (dist.f[dM00])[kN_M00];
+   fN[d0P0] = (dist.f[d0P0])[kN_000];
+   fN[d0M0] = (dist.f[d0M0])[kN_0M0];
+   fN[d00P] = (dist.f[d00P])[kN_000];
+   fN[d00M] = (dist.f[d00M])[kN_00M];
+   fN[dPP0] = (dist.f[dPP0])[kN_000];
+   fN[dMM0] = (dist.f[dMM0])[kN_MM0];
+   fN[dPM0] = (dist.f[dPM0])[kN_0M0];
+   fN[dMP0] = (dist.f[dMP0])[kN_M00];
+   fN[dP0P] = (dist.f[dP0P])[kN_000];
+   fN[dM0M] = (dist.f[dM0M])[kN_M0M];
+   fN[dP0M] = (dist.f[dP0M])[kN_00M];
+   fN[dM0P] = (dist.f[dM0P])[kN_M00];
+   fN[d0PP] = (dist.f[d0PP])[kN_000];
+   fN[d0MM] = (dist.f[d0MM])[kN_0MM];
+   fN[d0PM] = (dist.f[d0PM])[kN_00M];
+   fN[d0MP] = (dist.f[d0MP])[kN_0M0];
+   fN[dPPP] = (dist.f[dPPP])[kN_000];
+   fN[dMPP] = (dist.f[dMPP])[kN_M00];
+   fN[dPMP] = (dist.f[dPMP])[kN_0M0];
+   fN[dMMP] = (dist.f[dMMP])[kN_MM0];
+   fN[dPPM] = (dist.f[dPPM])[kN_00M];
+   fN[dMPM] = (dist.f[dMPM])[kN_M0M];
+   fN[dPMM] = (dist.f[dPMM])[kN_0MM];
+   fN[dMMM] = (dist.f[dMMM])[kN_MMM];
+   //////////////////////////////////////////////////////////////////////////
+   real drho = vf::lbm::getDensity(f);
+
+   real rhoCorrection = densityCorrectionFactor*drho;
+
+   real cs = c1o1 / sqrtf(c3o1);
+
+   getPointersToDistributions(dist, distributions, numberOfLBnodes, !isEvenTimestep);
+
+   switch(direction)
+   {
+      case dM00:
+         (dist.f[dP00])[k_000] = computeOutflowDistribution(f, fN, dP00  , rhoCorrection, cs, c2o27);
+         (dist.f[dPM0])[k_0M0] = computeOutflowDistribution(f, fN, dPM0, rhoCorrection, cs, c1o54);
+         (dist.f[dPP0])[k_000] = computeOutflowDistribution(f, fN, dPP0, rhoCorrection, cs, c1o54);
+         (dist.f[dP0M])[k_00M] = computeOutflowDistribution(f, fN, dP0M, rhoCorrection, cs, c1o54);
+         (dist.f[dP0P])[k_000] = computeOutflowDistribution(f, fN, dP0P, rhoCorrection, cs, c1o54);
+         (dist.f[dPMP])[k_0M0] = computeOutflowDistribution(f, fN, dPMP, rhoCorrection, cs, c1o216);
+         (dist.f[dPPP])[k_000] = computeOutflowDistribution(f, fN, dPPP, rhoCorrection, cs, c1o216);
+         (dist.f[dPMM])[k_0MM] = computeOutflowDistribution(f, fN, dPMM, rhoCorrection, cs, c1o216);
+         (dist.f[dPPM])[k_00M] = computeOutflowDistribution(f, fN, dPPM, rhoCorrection, cs, c1o216);
+         break;
+
+      case dP00:
+         (dist.f[dM00])[k_M00] = computeOutflowDistribution(f, fN, dM00, rhoCorrection, cs, c2o27);
+         (dist.f[dMM0])[k_MM0] = computeOutflowDistribution(f, fN, dMM0, rhoCorrection, cs, c1o54);
+         (dist.f[dMP0])[k_M00] = computeOutflowDistribution(f, fN, dMP0, rhoCorrection, cs, c1o54);
+         (dist.f[dM0M])[k_M0M] = computeOutflowDistribution(f, fN, dM0M, rhoCorrection, cs, c1o54);
+         (dist.f[dM0P])[k_M00] = computeOutflowDistribution(f, fN, dM0P, rhoCorrection, cs, c1o54);
+         (dist.f[dMMP])[k_MM0] = computeOutflowDistribution(f, fN, dMMP, rhoCorrection, cs, c1o216);
+         (dist.f[dMPP])[k_M00] = computeOutflowDistribution(f, fN, dMPP, rhoCorrection, cs, c1o216);
+         (dist.f[dMMM])[k_MMM] = computeOutflowDistribution(f, fN, dMMM, rhoCorrection, cs, c1o216);
+         (dist.f[dMPM])[k_M0M] = computeOutflowDistribution(f, fN, dMPM, rhoCorrection, cs, c1o216);
+         break;
+
+      case d0M0:
+         (dist.f[d0P0])[k_000] = computeOutflowDistribution(f, fN, d0P0, rhoCorrection, cs, c2o27);
+         (dist.f[dPP0])[k_000] = computeOutflowDistribution(f, fN, dPP0, rhoCorrection, cs, c1o54);
+         (dist.f[dMP0])[k_M00] = computeOutflowDistribution(f, fN, dMP0, rhoCorrection, cs, c1o54);
+         (dist.f[d0PP])[k_000] = computeOutflowDistribution(f, fN, d0PP, rhoCorrection, cs, c1o54);
+         (dist.f[d0PM])[k_00M] = computeOutflowDistribution(f, fN, d0PM, rhoCorrection, cs, c1o54);
+         (dist.f[dPPP])[k_000] = computeOutflowDistribution(f, fN, dPPP, rhoCorrection, cs, c1o216);
+         (dist.f[dMPP])[k_M00] = computeOutflowDistribution(f, fN, dMPP, rhoCorrection, cs, c1o216);
+         (dist.f[dPPM])[k_00M] = computeOutflowDistribution(f, fN, dPPM, rhoCorrection, cs, c1o216);
+         (dist.f[dMPM])[k_M0M] = computeOutflowDistribution(f, fN, dMPM, rhoCorrection, cs, c1o216);
+         break;
+
+      case d0P0:
+         (dist.f[d0M0])[k_0M0] =computeOutflowDistribution(f, fN, d0M0, rhoCorrection, cs, c2o27);
+         (dist.f[dPM0])[k_0M0] =computeOutflowDistribution(f, fN, dPM0, rhoCorrection, cs, c1o54);
+         (dist.f[dMM0])[k_MM0] =computeOutflowDistribution(f, fN, dMM0, rhoCorrection, cs, c1o54);
+         (dist.f[d0MP])[k_0M0] =computeOutflowDistribution(f, fN, d0MP, rhoCorrection, cs, c1o54);
+         (dist.f[d0MM])[k_0MM] =computeOutflowDistribution(f, fN, d0MM, rhoCorrection, cs, c1o54);
+         (dist.f[dPMP])[k_0M0] =computeOutflowDistribution(f, fN, dPMP, rhoCorrection, cs, c1o216);
+         (dist.f[dMMP])[k_MM0] =computeOutflowDistribution(f, fN, dMMP, rhoCorrection, cs, c1o216);
+         (dist.f[dPMM])[k_0MM] =computeOutflowDistribution(f, fN, dPMM, rhoCorrection, cs, c1o216);
+         (dist.f[dMMM])[k_MMM] =computeOutflowDistribution(f, fN, dMMM, rhoCorrection, cs, c1o216);
+         break;
+
+      case d00M:
+         (dist.f[d00P])[k_000] = computeOutflowDistribution(f, fN, d00P, rhoCorrection, cs, c2o27);
+         (dist.f[dP0P])[k_000] = computeOutflowDistribution(f, fN, dP0P, rhoCorrection, cs, c1o54);
+         (dist.f[dM0P])[k_M00] = computeOutflowDistribution(f, fN, dM0P, rhoCorrection, cs, c1o54);
+         (dist.f[d0PP])[k_000] = computeOutflowDistribution(f, fN, d0PP, rhoCorrection, cs, c1o54);
+         (dist.f[d0MP])[k_0M0] = computeOutflowDistribution(f, fN, d0MP, rhoCorrection, cs, c1o54);
+         (dist.f[dPPP])[k_000] = computeOutflowDistribution(f, fN, dPPP, rhoCorrection, cs, c1o216);
+         (dist.f[dMPP])[k_M00] = computeOutflowDistribution(f, fN, dMPP, rhoCorrection, cs, c1o216);
+         (dist.f[dPMP])[k_0M0] = computeOutflowDistribution(f, fN, dPMP, rhoCorrection, cs, c1o216);
+         (dist.f[dMMP])[k_MM0] = computeOutflowDistribution(f, fN, dMMP, rhoCorrection, cs, c1o216);
+         break;
+
+      case d00P:
+         (dist.f[d00M])[k_00M] = computeOutflowDistribution(f, fN, d00M, rhoCorrection, cs, c2o27);
+         (dist.f[dP0M])[k_00M] = computeOutflowDistribution(f, fN, dP0M, rhoCorrection, cs, c1o54);
+         (dist.f[dM0M])[k_M0M] = computeOutflowDistribution(f, fN, dM0M, rhoCorrection, cs, c1o54);
+         (dist.f[d0PM])[k_00M] = computeOutflowDistribution(f, fN, d0PM, rhoCorrection, cs, c1o54);
+         (dist.f[d0MM])[k_0MM] = computeOutflowDistribution(f, fN, d0MM, rhoCorrection, cs, c1o54);
+         (dist.f[dPPM])[k_00M] = computeOutflowDistribution(f, fN, dPPM, rhoCorrection, cs, c1o216);
+         (dist.f[dMPM])[k_M0M] = computeOutflowDistribution(f, fN, dMPM, rhoCorrection, cs, c1o216);
+         (dist.f[dPMM])[k_0MM] = computeOutflowDistribution(f, fN, dPMM, rhoCorrection, cs, c1o216);
+         (dist.f[dMMM])[k_MMM] = computeOutflowDistribution(f, fN, dMMM, rhoCorrection, cs, c1o216);
+         break;
+      default:
+         break;
+   }
+}
+
diff --git a/src/gpu/core/BoundaryConditions/Outflow/Outflow_Device.cuh b/src/gpu/core/BoundaryConditions/Outflow/Outflow_Device.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..c1039d03b6df127ffe670de0c53d1a94a69d7524
--- /dev/null
+++ b/src/gpu/core/BoundaryConditions/Outflow/Outflow_Device.cuh
@@ -0,0 +1,46 @@
+//  _    ___      __              __________      _     __        ______________   __
+// | |  / (_)____/ /___  ______ _/ / ____/ /_  __(_)___/ /____   /  ___/ __  / /  / /
+// | | / / / ___/ __/ / / / __ `/ / /_  / / / / / / __  / ___/  / /___/ /_/ / /  / /
+// | |/ / / /  / /_/ /_/ / /_/ / / __/ / / /_/ / / /_/ (__  )  / /_) / ____/ /__/ /
+// |___/_/_/   \__/\__,_/\__,_/_/_/   /_/\__,_/_/\__,_/____/   \____/_/    \_____/
+//
+//////////////////////////////////////////////////////////////////////////
+#ifndef Outflow_Device_H
+#define Outflow_Device_H
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+
+#include "LBM/LB.h"
+
+__global__ void OutflowNonReflecting_Device(  
+    real* rhoBC,
+    real* distributions,
+    int* k_Q,
+    int* k_N,
+    int numberOfBCnodes,
+    real om1,
+    unsigned int* neighborX,
+    unsigned int* neighborY,
+    unsigned int* neighborZ,
+    unsigned long long numberOfLBnodes,
+    bool isEvenTimestep,
+    int direction);
+
+__global__ void OutflowNonReflectingPressureCorrection_Device(  
+    real* rhoBC,
+    real* distributions,
+    int* k_Q,
+    int* k_N,
+    int numberOfBCnodes,
+    real om1,
+    unsigned int* neighborX,
+    unsigned int* neighborY,
+    unsigned int* neighborZ,
+    unsigned long long numberOfLBnodes,
+    bool isEvenTimestep,
+    int direction,
+    real densityCorrectionFactor);
+
+
+#endif
diff --git a/src/gpu/core/GPU/GPU_Interface.h b/src/gpu/core/GPU/GPU_Interface.h
index 2f25f5f15cb5f9334c1cef070ba9c32a7c11c68b..6d74a76f0f5b62ce6b9cb896478b394e37b8c662 100644
--- a/src/gpu/core/GPU/GPU_Interface.h
+++ b/src/gpu/core/GPU/GPU_Interface.h
@@ -576,10 +576,6 @@ void QPressDevDirDepBot27(unsigned int numberOfThreads,
                                      unsigned long long numberOfLBnodes, 
                                      bool isEvenTimestep);
 
-void QPressNoRhoDev27(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition);
-
-void QPressZeroRhoOutflowDev27(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition);
-
 void QInflowScaleByPressDev27(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition);
 
 void QPressDevOld27(unsigned int numberOfThreads,
diff --git a/src/gpu/core/GPU/GPU_Kernels.cuh b/src/gpu/core/GPU/GPU_Kernels.cuh
index 6e007fb6ced0b4301b0f212d4b0199ae4988b728..4c554d6efba903ff7d3710f5b351897908e6420f 100644
--- a/src/gpu/core/GPU/GPU_Kernels.cuh
+++ b/src/gpu/core/GPU/GPU_Kernels.cuh
@@ -781,33 +781,6 @@ __global__ void QPressDeviceDirDepBot27(  real* rhoBC,
                                                      unsigned long long numberOfLBnodes,
                                                      bool isEvenTimestep);
 
-__global__ void QPressNoRhoDevice27(  real* rhoBC,
-                                                 real* distributions,
-                                                 int* k_Q,
-                                                 int* k_N,
-                                                 int numberOfBCnodes,
-                                                 real om1,
-                                                 unsigned int* neighborX,
-                                                 unsigned int* neighborY,
-                                                 unsigned int* neighborZ,
-                                                 unsigned long long numberOfLBnodes,
-                                                 bool isEvenTimestep,
-                                                 int direction);
-
-__global__ void QPressZeroRhoOutflowDevice27(  real* rhoBC,
-                                            real* distributions,
-                                            int* k_Q,
-                                            int* k_N,
-                                            int numberOfBCnodes,
-                                            real om1,
-                                            unsigned int* neighborX,
-                                            unsigned int* neighborY,
-                                            unsigned int* neighborZ,
-                                            unsigned long long numberOfLBnodes,
-                                            bool isEvenTimestep,
-                                            int direction,
-                                            real densityCorrectionFactor);
-
 __global__ void QInflowScaleByPressDevice27(  real* rhoBC,
                                                          real* DD,
                                                          int* k_Q,
diff --git a/src/gpu/core/GPU/LBMKernel.cu b/src/gpu/core/GPU/LBMKernel.cu
index a870d43bb36811441afa654a01f4eeb7671ba28e..d5837f86340c45f79ca8cedb57da11d20ae63842 100644
--- a/src/gpu/core/GPU/LBMKernel.cu
+++ b/src/gpu/core/GPU/LBMKernel.cu
@@ -2208,49 +2208,6 @@ void QPressDevDirDepBot27(
     getLastCudaError("QPressDeviceDirDepBot27 execution failed");
 }
 //////////////////////////////////////////////////////////////////////////
-void QPressNoRhoDev27(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition)
-{
-    dim3 grid = vf::cuda::getCudaGrid( parameterDevice->numberofthreads,  boundaryCondition->numberOfBCnodes);
-    dim3 threads(parameterDevice->numberofthreads, 1, 1 );
-
-    QPressNoRhoDevice27<<< grid, threads >>> (
-        boundaryCondition->RhoBC,
-        parameterDevice->distributions.f[0],
-        boundaryCondition->k,
-        boundaryCondition->kN,
-        boundaryCondition->numberOfBCnodes,
-        parameterDevice->omega,
-        parameterDevice->neighborX,
-        parameterDevice->neighborY,
-        parameterDevice->neighborZ,
-        parameterDevice->numberOfNodes,
-        parameterDevice->isEvenTimestep,
-        vf::lbm::dir::dP00);
-    getLastCudaError("QPressNoRhoDevice27 execution failed");
-}
-//////////////////////////////////////////////////////////////////////////
-void QPressZeroRhoOutflowDev27(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition)
-{
-    dim3 grid = vf::cuda::getCudaGrid( parameterDevice->numberofthreads,  boundaryCondition->numberOfBCnodes);
-    dim3 threads(parameterDevice->numberofthreads, 1, 1 );
-
-    QPressZeroRhoOutflowDevice27<<< grid, threads >>> (
-        boundaryCondition->RhoBC,
-        parameterDevice->distributions.f[0],
-        boundaryCondition->k,
-        boundaryCondition->kN,
-        boundaryCondition->numberOfBCnodes,
-        parameterDevice->omega,
-        parameterDevice->neighborX,
-        parameterDevice->neighborY,
-        parameterDevice->neighborZ,
-        parameterDevice->numberOfNodes,
-        parameterDevice->isEvenTimestep,
-        vf::lbm::dir::dP00,
-        parameterDevice->outflowPressureCorrectionFactor);
-    getLastCudaError("QPressZeroRhoOutflowDevice27 execution failed");
-}
-//////////////////////////////////////////////////////////////////////////
 void QInflowScaleByPressDev27(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition)
 {
     dim3 grid = vf::cuda::getCudaGrid( parameterDevice->numberofthreads,  boundaryCondition->numberOfBCnodes);
diff --git a/src/gpu/core/GPU/PressBCs27.cu b/src/gpu/core/GPU/PressBCs27.cu
index 2c2a4f0380d8d8b7e6856292894d298c228ada28..a146e34a5ca681d3359d84c4cbd77dae84729cea 100644
--- a/src/gpu/core/GPU/PressBCs27.cu
+++ b/src/gpu/core/GPU/PressBCs27.cu
@@ -2812,508 +2812,6 @@ __global__ void QPressDeviceDirDepBot27(
 
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-__host__ __device__ real computeOutflowDistribution(const real* const &f, const real* const &f1, const int dir, const real cs)
-{
-   return f1[dir] * cs + (c1o1 - cs) * f[dir];
-}
-
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-__global__ void QPressNoRhoDevice27(
-    real* rhoBC,
-    real* distributions,
-    int* k_Q,
-    int* k_N,
-    int numberOfBCnodes,
-    real om1,
-    unsigned int* neighborX,
-    unsigned int* neighborY,
-    unsigned int* neighborZ,
-    unsigned long long numberOfLBnodes,
-    bool isEvenTimestep,
-    int direction)
-{
-   ////////////////////////////////////////////////////////////////////////////////
-   //! - Get the node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
-   //!
-   const unsigned nodeIndex = getNodeIndex();
-
-   //////////////////////////////////////////////////////////////////////////
-
-   if(nodeIndex >= numberOfBCnodes) return;
-
-   ////////////////////////////////////////////////////////////////////////////////
-   //index
-   unsigned int KQK  = k_Q[nodeIndex];
-   // unsigned int kzero= KQK;
-   unsigned int ke   = KQK;
-   unsigned int kw   = neighborX[KQK];
-   unsigned int kn   = KQK;
-   unsigned int ks   = neighborY[KQK];
-   unsigned int kt   = KQK;
-   unsigned int kb   = neighborZ[KQK];
-   unsigned int ksw  = neighborY[kw];
-   unsigned int kne  = KQK;
-   unsigned int kse  = ks;
-   unsigned int knw  = kw;
-   unsigned int kbw  = neighborZ[kw];
-   unsigned int kte  = KQK;
-   unsigned int kbe  = kb;
-   unsigned int ktw  = kw;
-   unsigned int kbs  = neighborZ[ks];
-   unsigned int ktn  = KQK;
-   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 = KQK;
-   unsigned int kbsw = neighborZ[ksw];
-   ////////////////////////////////////////////////////////////////////////////////
-   //index1
-   unsigned int K1QK  = k_N[nodeIndex];
-   //unsigned int k1zero= K1QK;
-   unsigned int k1e   = K1QK;
-   unsigned int k1w   = neighborX[K1QK];
-   unsigned int k1n   = K1QK;
-   unsigned int k1s   = neighborY[K1QK];
-   unsigned int k1t   = K1QK;
-   unsigned int k1b   = neighborZ[K1QK];
-   unsigned int k1sw  = neighborY[k1w];
-   unsigned int k1ne  = K1QK;
-   unsigned int k1se  = k1s;
-   unsigned int k1nw  = k1w;
-   unsigned int k1bw  = neighborZ[k1w];
-   unsigned int k1te  = K1QK;
-   unsigned int k1be  = k1b;
-   unsigned int k1tw  = k1w;
-   unsigned int k1bs  = neighborZ[k1s];
-   unsigned int k1tn  = K1QK;
-   unsigned int k1bn  = k1b;
-   unsigned int k1ts  = k1s;
-   unsigned int k1tse = k1s;
-   unsigned int k1bnw = k1bw;
-   unsigned int k1tnw = k1w;
-   unsigned int k1bse = k1bs;
-   unsigned int k1tsw = k1sw;
-   unsigned int k1bne = k1b;
-   unsigned int k1tne = K1QK;
-   unsigned int k1bsw = neighborZ[k1sw];
-   ////////////////////////////////////////////////////////////////////////////////
-   Distributions27 dist;
-   getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
-   real f[27], f1[27];
-   //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-   f1[dP00] = (dist.f[dP00])[k1e   ];
-   f1[dM00] = (dist.f[dM00])[k1w   ];
-   f1[d0P0] = (dist.f[d0P0])[k1n   ];
-   f1[d0M0] = (dist.f[d0M0])[k1s   ];
-   f1[d00P] = (dist.f[d00P])[k1t   ];
-   f1[d00M] = (dist.f[d00M])[k1b   ];
-   f1[dPP0] = (dist.f[dPP0])[k1ne  ];
-   f1[dMM0] = (dist.f[dMM0])[k1sw  ];
-   f1[dPM0] = (dist.f[dPM0])[k1se  ];
-   f1[dMP0] = (dist.f[dMP0])[k1nw  ];
-   f1[dP0P] = (dist.f[dP0P])[k1te  ];
-   f1[dM0M] = (dist.f[dM0M])[k1bw  ];
-   f1[dP0M] = (dist.f[dP0M])[k1be  ];
-   f1[dM0P] = (dist.f[dM0P])[k1tw  ];
-   f1[d0PP] = (dist.f[d0PP])[k1tn  ];
-   f1[d0MM] = (dist.f[d0MM])[k1bs  ];
-   f1[d0PM] = (dist.f[d0PM])[k1bn  ];
-   f1[d0MP] = (dist.f[d0MP])[k1ts  ];
-   // f1[d000] = (dist.f[d000])[k1zero];
-   f1[dPPP] = (dist.f[dPPP])[k1tne ];
-   f1[dMMP] = (dist.f[dMMP])[k1tsw ];
-   f1[dPMP] = (dist.f[dPMP])[k1tse ];
-   f1[dMPP] = (dist.f[dMPP])[k1tnw ];
-   f1[dPPM] = (dist.f[dPPM])[k1bne ];
-   f1[dMMM] = (dist.f[dMMM])[k1bsw ];
-   f1[dPMM] = (dist.f[dPMM])[k1bse ];
-   f1[dMPM] = (dist.f[dMPM])[k1bnw ];
-   //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-   f[dP00] = (dist.f[dP00])[ke   ];
-   f[dM00] = (dist.f[dM00])[kw   ];
-   f[d0P0] = (dist.f[d0P0])[kn   ];
-   f[d0M0] = (dist.f[d0M0])[ks   ];
-   f[d00P] = (dist.f[d00P])[kt   ];
-   f[d00M] = (dist.f[d00M])[kb   ];
-   f[dPP0] = (dist.f[dPP0])[kne  ];
-   f[dMM0] = (dist.f[dMM0])[ksw  ];
-   f[dPM0] = (dist.f[dPM0])[kse  ];
-   f[dMP0] = (dist.f[dMP0])[knw  ];
-   f[dP0P] = (dist.f[dP0P])[kte  ];
-   f[dM0M] = (dist.f[dM0M])[kbw  ];
-   f[dP0M] = (dist.f[dP0M])[kbe  ];
-   f[dM0P] = (dist.f[dM0P])[ktw  ];
-   f[d0PP] = (dist.f[d0PP])[ktn  ];
-   f[d0MM] = (dist.f[d0MM])[kbs  ];
-   f[d0PM] = (dist.f[d0PM])[kbn  ];
-   f[d0MP] = (dist.f[d0MP])[kts  ];
-   // f[d000] = (dist.f[d000])[kzero];
-   f[dPPP] = (dist.f[dPPP])[ktne ];
-   f[dMMP] = (dist.f[dMMP])[ktsw ];
-   f[dPMP] = (dist.f[dPMP])[ktse ];
-   f[dMPP] = (dist.f[dMPP])[ktnw ];
-   f[dPPM] = (dist.f[dPPM])[kbne ];
-   f[dMMM] = (dist.f[dMMM])[kbsw ];
-   f[dPMM] = (dist.f[dPMM])[kbse ];
-   f[dMPM] = (dist.f[dMPM])[kbnw ];
-   //////////////////////////////////////////////////////////////////////////
-
-
-   real cs = c1o1 / sqrtf(c3o1);
-
-   //////////////////////////////////////////////////////////////////////////
-   getPointersToDistributions(dist, distributions, numberOfLBnodes, !isEvenTimestep);
-   switch(direction)
-   {
-      case dM00:
-         (dist.f[dP00])[ke   ] = computeOutflowDistribution(f, f1, dP00, cs);
-         (dist.f[dPM0])[kse  ] = computeOutflowDistribution(f, f1, dPM0, cs);
-         (dist.f[dPP0])[kne  ] = computeOutflowDistribution(f, f1, dPP0, cs);
-         (dist.f[dP0M])[kbe  ] = computeOutflowDistribution(f, f1, dP0M, cs);
-         (dist.f[dP0P])[kte  ] = computeOutflowDistribution(f, f1, dP0P, cs);
-         (dist.f[dPMP])[ktse ] = computeOutflowDistribution(f, f1, dPMP, cs);
-         (dist.f[dPPP])[ktne ] = computeOutflowDistribution(f, f1, dPPP, cs);
-         (dist.f[dPMM])[kbse ] = computeOutflowDistribution(f, f1, dPMM, cs);
-         (dist.f[dPPM])[kbne ] = computeOutflowDistribution(f, f1, dPPM, cs);
-         break;
-
-      case dP00:
-         (dist.f[dM00])[kw   ] = computeOutflowDistribution(f, f1, dM00, cs);
-         (dist.f[dMM0])[ksw  ] = computeOutflowDistribution(f, f1, dMM0, cs);
-         (dist.f[dMP0])[knw  ] = computeOutflowDistribution(f, f1, dMP0, cs);
-         (dist.f[dM0M])[kbw  ] = computeOutflowDistribution(f, f1, dM0M, cs);
-         (dist.f[dM0P])[ktw  ] = computeOutflowDistribution(f, f1, dM0P, cs);
-         (dist.f[dMMP])[ktsw ] = computeOutflowDistribution(f, f1, dMMP, cs);
-         (dist.f[dMPP])[ktnw ] = computeOutflowDistribution(f, f1, dMPP, cs);
-         (dist.f[dMMM])[kbsw ] = computeOutflowDistribution(f, f1, dMMM, cs);
-         (dist.f[dMPM])[kbnw ] = computeOutflowDistribution(f, f1, dMPM, cs);
-         break;
-
-      case d0M0:
-         (dist.f[d0P0])[kn   ] = computeOutflowDistribution(f, f1, d0P0, cs);
-         (dist.f[dPP0])[kne  ] = computeOutflowDistribution(f, f1, dPP0, cs);
-         (dist.f[dMP0])[knw  ] = computeOutflowDistribution(f, f1, dMP0, cs);
-         (dist.f[d0PP])[ktn  ] = computeOutflowDistribution(f, f1, d0PP, cs);
-         (dist.f[d0PM])[kbn  ] = computeOutflowDistribution(f, f1, d0PM, cs);
-         (dist.f[dPPP])[ktne ] = computeOutflowDistribution(f, f1, dPPP, cs);
-         (dist.f[dMPP])[ktnw ] = computeOutflowDistribution(f, f1, dMPP, cs);
-         (dist.f[dPPM])[kbne ] = computeOutflowDistribution(f, f1, dPPM, cs);
-         (dist.f[dMPM])[kbnw ] = computeOutflowDistribution(f, f1, dMPM, cs);
-         break;
-
-      case d0P0:
-         (dist.f[d0M0])[ks   ] = computeOutflowDistribution(f, f1, d0M0, cs);
-         (dist.f[dPM0])[kse  ] = computeOutflowDistribution(f, f1, dPM0, cs);
-         (dist.f[dMM0])[ksw  ] = computeOutflowDistribution(f, f1, dMM0, cs);
-         (dist.f[d0MP])[kts  ] = computeOutflowDistribution(f, f1, d0MP, cs);
-         (dist.f[d0MM])[kbs  ] = computeOutflowDistribution(f, f1, d0MM, cs);
-         (dist.f[dPMP])[ktse ] = computeOutflowDistribution(f, f1, dPMP, cs);
-         (dist.f[dMMP])[ktsw ] = computeOutflowDistribution(f, f1, dMMP, cs);
-         (dist.f[dPMM])[kbse ] = computeOutflowDistribution(f, f1, dPMM, cs);
-         (dist.f[dMMM])[kbsw ] = computeOutflowDistribution(f, f1, dMMM, cs);
-         break;
-
-      case d00M:
-         (dist.f[d00P])[kt   ] = computeOutflowDistribution(f, f1, d00P, cs);
-         (dist.f[dP0P])[kte  ] = computeOutflowDistribution(f, f1, dP0P, cs);
-         (dist.f[dM0P])[ktw  ] = computeOutflowDistribution(f, f1, dM0P, cs);
-         (dist.f[d0PP])[ktn  ] = computeOutflowDistribution(f, f1, d0PP, cs);
-         (dist.f[d0MP])[kts  ] = computeOutflowDistribution(f, f1, d0MP, cs);
-         (dist.f[dPPP])[ktne ] = computeOutflowDistribution(f, f1, dPPP, cs);
-         (dist.f[dMPP])[ktnw ] = computeOutflowDistribution(f, f1, dMPP, cs);
-         (dist.f[dPMP])[ktse ] = computeOutflowDistribution(f, f1, dPMP, cs);
-         (dist.f[dMMP])[ktsw ] = computeOutflowDistribution(f, f1, dMMP, cs);
-         break;
-
-      case d00P:
-         (dist.f[d00M])[kb   ] = computeOutflowDistribution(f, f1, d00M, cs);
-         (dist.f[dP0M])[kbe  ] = computeOutflowDistribution(f, f1, dP0M, cs);
-         (dist.f[dM0M])[kbw  ] = computeOutflowDistribution(f, f1, dM0M, cs);
-         (dist.f[d0PM])[kbn  ] = computeOutflowDistribution(f, f1, d0PM, cs);
-         (dist.f[d0MM])[kbs  ] = computeOutflowDistribution(f, f1, d0MM, cs);
-         (dist.f[dPPM])[kbne ] = computeOutflowDistribution(f, f1, dPPM, cs);
-         (dist.f[dMPM])[kbnw ] = computeOutflowDistribution(f, f1, dMPM, cs);
-         (dist.f[dPMM])[kbse ] = computeOutflowDistribution(f, f1, dPMM, cs);
-         (dist.f[dMMM])[kbsw ] = computeOutflowDistribution(f, f1, dMMM, cs);
-         break;
-      default:
-         break;
-   }
-}
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-__host__ __device__ real computeOutflowDistribution(const real* const &f, const real* const &f1, const int dir, const real rhoCorrection, const real cs, const real weight)
-{
-   return f1[dir  ] * cs + (c1o1 - cs) * f[dir  ] - weight *rhoCorrection;
-}
-
-__global__ void QPressZeroRhoOutflowDevice27(
-    real* rhoBC,
-    real* distributions,
-    int* k_Q,
-    int* k_N,
-    int numberOfBCnodes,
-    real om1,
-    unsigned int* neighborX,
-    unsigned int* neighborY,
-    unsigned int* neighborZ,
-    unsigned long long numberOfLBnodes,
-    bool isEvenTimestep,
-    int direction,
-    real densityCorrectionFactor)
-{
-   ////////////////////////////////////////////////////////////////////////////////
-   //! - Get the node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
-   //!
-   const unsigned nodeIndex = getNodeIndex();
-
-   //////////////////////////////////////////////////////////////////////////
-
-   if( nodeIndex >= numberOfBCnodes ) return;
-
-   ////////////////////////////////////////////////////////////////////////////////
-   //index
-
-   uint k_000 = k_Q[nodeIndex];
-   uint k_M00 = neighborX[k_000];
-   uint k_0M0 = neighborY[k_000];
-   uint k_00M = neighborZ[k_000];
-   uint k_MM0 = neighborY[k_M00];
-   uint k_M0M = neighborZ[k_M00];
-   uint k_0MM = neighborZ[k_0M0];
-   uint k_MMM = neighborZ[k_MM0];
-
-   ////////////////////////////////////////////////////////////////////////////////
-   //index of neighbor
-   uint kN_000 = k_N[nodeIndex];
-   uint kN_M00 = neighborX[k_000];
-   uint kN_0M0 = neighborY[k_000];
-   uint kN_00M = neighborZ[k_000];
-   uint kN_MM0 = neighborY[k_M00];
-   uint kN_M0M = neighborZ[k_M00];
-   uint kN_0MM = neighborZ[k_0M0];
-   uint kN_MMM = neighborZ[k_MM0];
-   ////////////////////////////////////////////////////////////////////////////////
-   Distributions27 dist;
-   getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
-   real f[27], fN[27];
-   //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-   f[d000] = (dist.f[d000])[k_000];
-   f[dP00] = (dist.f[dP00])[k_000];
-   f[dM00] = (dist.f[dM00])[k_M00];
-   f[d0P0] = (dist.f[d0P0])[k_000];
-   f[d0M0] = (dist.f[d0M0])[k_0M0];
-   f[d00P] = (dist.f[d00P])[k_000];
-   f[d00M] = (dist.f[d00M])[k_00M];
-   f[dPP0] = (dist.f[dPP0])[k_000];
-   f[dMM0] = (dist.f[dMM0])[k_MM0];
-   f[dPM0] = (dist.f[dPM0])[k_0M0];
-   f[dMP0] = (dist.f[dMP0])[k_M00];
-   f[dP0P] = (dist.f[dP0P])[k_000];
-   f[dM0M] = (dist.f[dM0M])[k_M0M];
-   f[dP0M] = (dist.f[dP0M])[k_00M];
-   f[dM0P] = (dist.f[dM0P])[k_M00];
-   f[d0PP] = (dist.f[d0PP])[k_000];
-   f[d0MM] = (dist.f[d0MM])[k_0MM];
-   f[d0PM] = (dist.f[d0PM])[k_00M];
-   f[d0MP] = (dist.f[d0MP])[k_0M0];
-   f[dPPP] = (dist.f[dPPP])[k_000];
-   f[dMPP] = (dist.f[dMPP])[k_M00];
-   f[dPMP] = (dist.f[dPMP])[k_0M0];
-   f[dMMP] = (dist.f[dMMP])[k_MM0];
-   f[dPPM] = (dist.f[dPPM])[k_00M];
-   f[dMPM] = (dist.f[dMPM])[k_M0M];
-   f[dPMM] = (dist.f[dPMM])[k_0MM];
-   f[dMMM] = (dist.f[dMMM])[k_MMM];
-   //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-   fN[d000] = (dist.f[d000])[kN_000];
-   fN[dP00] = (dist.f[dP00])[kN_000];
-   fN[dM00] = (dist.f[dM00])[kN_M00];
-   fN[d0P0] = (dist.f[d0P0])[kN_000];
-   fN[d0M0] = (dist.f[d0M0])[kN_0M0];
-   fN[d00P] = (dist.f[d00P])[kN_000];
-   fN[d00M] = (dist.f[d00M])[kN_00M];
-   fN[dPP0] = (dist.f[dPP0])[kN_000];
-   fN[dMM0] = (dist.f[dMM0])[kN_MM0];
-   fN[dPM0] = (dist.f[dPM0])[kN_0M0];
-   fN[dMP0] = (dist.f[dMP0])[kN_M00];
-   fN[dP0P] = (dist.f[dP0P])[kN_000];
-   fN[dM0M] = (dist.f[dM0M])[kN_M0M];
-   fN[dP0M] = (dist.f[dP0M])[kN_00M];
-   fN[dM0P] = (dist.f[dM0P])[kN_M00];
-   fN[d0PP] = (dist.f[d0PP])[kN_000];
-   fN[d0MM] = (dist.f[d0MM])[kN_0MM];
-   fN[d0PM] = (dist.f[d0PM])[kN_00M];
-   fN[d0MP] = (dist.f[d0MP])[kN_0M0];
-   fN[dPPP] = (dist.f[dPPP])[kN_000];
-   fN[dMPP] = (dist.f[dMPP])[kN_M00];
-   fN[dPMP] = (dist.f[dPMP])[kN_0M0];
-   fN[dMMP] = (dist.f[dMMP])[kN_MM0];
-   fN[dPPM] = (dist.f[dPPM])[kN_00M];
-   fN[dMPM] = (dist.f[dMPM])[kN_M0M];
-   fN[dPMM] = (dist.f[dPMM])[kN_0MM];
-   fN[dMMM] = (dist.f[dMMM])[kN_MMM];
-   //////////////////////////////////////////////////////////////////////////
-   real drho = vf::lbm::getDensity(f);
-
-   real rhoCorrection = densityCorrectionFactor*drho;
-
-   real cs = c1o1 / sqrtf(c3o1);
-
-   getPointersToDistributions(dist, distributions, numberOfLBnodes, !isEvenTimestep);
-
-   switch(direction)
-   {
-      case dM00:
-         (dist.f[dP00])[k_000] = computeOutflowDistribution(f, fN, dP00  , rhoCorrection, cs, c2o27);
-         (dist.f[dPM0])[k_0M0] = computeOutflowDistribution(f, fN, dPM0, rhoCorrection, cs, c1o54);
-         (dist.f[dPP0])[k_000] = computeOutflowDistribution(f, fN, dPP0, rhoCorrection, cs, c1o54);
-         (dist.f[dP0M])[k_00M] = computeOutflowDistribution(f, fN, dP0M, rhoCorrection, cs, c1o54);
-         (dist.f[dP0P])[k_000] = computeOutflowDistribution(f, fN, dP0P, rhoCorrection, cs, c1o54);
-         (dist.f[dPMP])[k_0M0] = computeOutflowDistribution(f, fN, dPMP, rhoCorrection, cs, c1o216);
-         (dist.f[dPPP])[k_000] = computeOutflowDistribution(f, fN, dPPP, rhoCorrection, cs, c1o216);
-         (dist.f[dPMM])[k_0MM] = computeOutflowDistribution(f, fN, dPMM, rhoCorrection, cs, c1o216);
-         (dist.f[dPPM])[k_00M] = computeOutflowDistribution(f, fN, dPPM, rhoCorrection, cs, c1o216);
-         break;
-
-      case dP00:
-         (dist.f[dM00])[k_M00] = computeOutflowDistribution(f, fN, dM00, rhoCorrection, cs, c2o27);
-         (dist.f[dMM0])[k_MM0] = computeOutflowDistribution(f, fN, dMM0, rhoCorrection, cs, c1o54);
-         (dist.f[dMP0])[k_M00] = computeOutflowDistribution(f, fN, dMP0, rhoCorrection, cs, c1o54);
-         (dist.f[dM0M])[k_M0M] = computeOutflowDistribution(f, fN, dM0M, rhoCorrection, cs, c1o54);
-         (dist.f[dM0P])[k_M00] = computeOutflowDistribution(f, fN, dM0P, rhoCorrection, cs, c1o54);
-         (dist.f[dMMP])[k_MM0] = computeOutflowDistribution(f, fN, dMMP, rhoCorrection, cs, c1o216);
-         (dist.f[dMPP])[k_M00] = computeOutflowDistribution(f, fN, dMPP, rhoCorrection, cs, c1o216);
-         (dist.f[dMMM])[k_MMM] = computeOutflowDistribution(f, fN, dMMM, rhoCorrection, cs, c1o216);
-         (dist.f[dMPM])[k_M0M] = computeOutflowDistribution(f, fN, dMPM, rhoCorrection, cs, c1o216);
-         break;
-
-      case d0M0:
-         (dist.f[d0P0])[k_000] = computeOutflowDistribution(f, fN, d0P0, rhoCorrection, cs, c2o27);
-         (dist.f[dPP0])[k_000] = computeOutflowDistribution(f, fN, dPP0, rhoCorrection, cs, c1o54);
-         (dist.f[dMP0])[k_M00] = computeOutflowDistribution(f, fN, dMP0, rhoCorrection, cs, c1o54);
-         (dist.f[d0PP])[k_000] = computeOutflowDistribution(f, fN, d0PP, rhoCorrection, cs, c1o54);
-         (dist.f[d0PM])[k_00M] = computeOutflowDistribution(f, fN, d0PM, rhoCorrection, cs, c1o54);
-         (dist.f[dPPP])[k_000] = computeOutflowDistribution(f, fN, dPPP, rhoCorrection, cs, c1o216);
-         (dist.f[dMPP])[k_M00] = computeOutflowDistribution(f, fN, dMPP, rhoCorrection, cs, c1o216);
-         (dist.f[dPPM])[k_00M] = computeOutflowDistribution(f, fN, dPPM, rhoCorrection, cs, c1o216);
-         (dist.f[dMPM])[k_M0M] = computeOutflowDistribution(f, fN, dMPM, rhoCorrection, cs, c1o216);
-         break;
-
-      case d0P0:
-         (dist.f[d0M0])[k_0M0] =computeOutflowDistribution(f, fN, d0M0, rhoCorrection, cs, c2o27);
-         (dist.f[dPM0])[k_0M0] =computeOutflowDistribution(f, fN, dPM0, rhoCorrection, cs, c1o54);
-         (dist.f[dMM0])[k_MM0] =computeOutflowDistribution(f, fN, dMM0, rhoCorrection, cs, c1o54);
-         (dist.f[d0MP])[k_0M0] =computeOutflowDistribution(f, fN, d0MP, rhoCorrection, cs, c1o54);
-         (dist.f[d0MM])[k_0MM] =computeOutflowDistribution(f, fN, d0MM, rhoCorrection, cs, c1o54);
-         (dist.f[dPMP])[k_0M0] =computeOutflowDistribution(f, fN, dPMP, rhoCorrection, cs, c1o216);
-         (dist.f[dMMP])[k_MM0] =computeOutflowDistribution(f, fN, dMMP, rhoCorrection, cs, c1o216);
-         (dist.f[dPMM])[k_0MM] =computeOutflowDistribution(f, fN, dPMM, rhoCorrection, cs, c1o216);
-         (dist.f[dMMM])[k_MMM] =computeOutflowDistribution(f, fN, dMMM, rhoCorrection, cs, c1o216);
-         break;
-
-      case d00M:
-         (dist.f[d00P])[k_000] = computeOutflowDistribution(f, fN, d00P, rhoCorrection, cs, c2o27);
-         (dist.f[dP0P])[k_000] = computeOutflowDistribution(f, fN, dP0P, rhoCorrection, cs, c1o54);
-         (dist.f[dM0P])[k_M00] = computeOutflowDistribution(f, fN, dM0P, rhoCorrection, cs, c1o54);
-         (dist.f[d0PP])[k_000] = computeOutflowDistribution(f, fN, d0PP, rhoCorrection, cs, c1o54);
-         (dist.f[d0MP])[k_0M0] = computeOutflowDistribution(f, fN, d0MP, rhoCorrection, cs, c1o54);
-         (dist.f[dPPP])[k_000] = computeOutflowDistribution(f, fN, dPPP, rhoCorrection, cs, c1o216);
-         (dist.f[dMPP])[k_M00] = computeOutflowDistribution(f, fN, dMPP, rhoCorrection, cs, c1o216);
-         (dist.f[dPMP])[k_0M0] = computeOutflowDistribution(f, fN, dPMP, rhoCorrection, cs, c1o216);
-         (dist.f[dMMP])[k_MM0] = computeOutflowDistribution(f, fN, dMMP, rhoCorrection, cs, c1o216);
-         break;
-
-      case d00P:
-         (dist.f[d00M])[k_00M] = computeOutflowDistribution(f, fN, d00M, rhoCorrection, cs, c2o27);
-         (dist.f[dP0M])[k_00M] = computeOutflowDistribution(f, fN, dP0M, rhoCorrection, cs, c1o54);
-         (dist.f[dM0M])[k_M0M] = computeOutflowDistribution(f, fN, dM0M, rhoCorrection, cs, c1o54);
-         (dist.f[d0PM])[k_00M] = computeOutflowDistribution(f, fN, d0PM, rhoCorrection, cs, c1o54);
-         (dist.f[d0MM])[k_0MM] = computeOutflowDistribution(f, fN, d0MM, rhoCorrection, cs, c1o54);
-         (dist.f[dPPM])[k_00M] = computeOutflowDistribution(f, fN, dPPM, rhoCorrection, cs, c1o216);
-         (dist.f[dMPM])[k_M0M] = computeOutflowDistribution(f, fN, dMPM, rhoCorrection, cs, c1o216);
-         (dist.f[dPMM])[k_0MM] = computeOutflowDistribution(f, fN, dPMM, rhoCorrection, cs, c1o216);
-         (dist.f[dMMM])[k_MMM] = computeOutflowDistribution(f, fN, dMMM, rhoCorrection, cs, c1o216);
-         break;
-      default:
-         break;
-   }
-}
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-