diff --git a/.gitignore b/.gitignore
index ae19800aa7dd8a859144426e8280eb3718add848..f87c8efbbd3b3877bd77212d6c2184db2aa409f1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -24,7 +24,7 @@ output/
 logs/
 
 # grid
-grid/
+.grid/
 
 # scripts
 scripts/
diff --git a/apps/gpu/HULC/main.cpp b/apps/gpu/HULC/main.cpp
index 35d70d1f49b344e0e5530c46582940ec581e3c7f..0c0be6088470ff9cd72baa14544209f01c024cb9 100644
--- a/apps/gpu/HULC/main.cpp
+++ b/apps/gpu/HULC/main.cpp
@@ -90,8 +90,8 @@ void setParameters(std::shared_ptr<Parameter> para, std::unique_ptr<input::Input
     para->setTemperatureInit(StringUtil::toFloat(input->getValue("Temp")));
     para->setTemperatureBC(StringUtil::toFloat(input->getValue("TempBC")));
     //////////////////////////////////////////////////////////////////////////
-    para->setViscosity(StringUtil::toFloat(input->getValue("Viscosity_LB")));
-    para->setVelocity(StringUtil::toFloat(input->getValue("Velocity_LB")));
+    para->setViscosityLB(StringUtil::toFloat(input->getValue("Viscosity_LB")));
+    para->setVelocityLB(StringUtil::toFloat(input->getValue("Velocity_LB")));
     para->setViscosityRatio(StringUtil::toFloat(input->getValue("Viscosity_Ratio_World_to_LB")));
     para->setVelocityRatio(StringUtil::toFloat(input->getValue("Velocity_Ratio_World_to_LB")));
     para->setDensityRatio(StringUtil::toFloat(input->getValue("Density_Ratio_World_to_LB")));
diff --git a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
index 9a44362d1b0753736e719aeefa100890174ad799..51ebaf7e465d75f793ba04d3ecd6686aa94b8b8d 100644
--- a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
+++ b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
@@ -145,8 +145,8 @@ void multipleLevel(const std::string& configPath)
     para->setMaxLevel(1);
 
 
-    para->setVelocity(velocityLB);
-    para->setViscosity(viscosityLB);
+    para->setVelocityLB(velocityLB);
+    para->setViscosityLB(viscosityLB);
     para->setVelocityRatio( dx / dt );
     para->setViscosityRatio( dx*dx/dt );
     para->setMainKernel("CumulantK17CompChim");
diff --git a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
index 781ac3be1f011c7c967821a7f5784e7020663c44..f6b2c248de8bcf7a76b4b8e6504f72d36a0d6eb8 100644
--- a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
+++ b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
@@ -157,8 +157,8 @@ void multipleLevel(const std::string& configPath)
     para->setPrintFiles(true);
 
     para->setForcing(pressureGradientLB, 0, 0);
-    para->setVelocity(velocityLB);
-    para->setViscosity(viscosityLB);
+    para->setVelocityLB(velocityLB);
+    para->setViscosityLB(viscosityLB);
     para->setVelocityRatio( dx / dt );
     para->setViscosityRatio( dx*dx/dt );
     para->setDensityRatio( 1.0 );
diff --git a/apps/gpu/LBM/DrivenCavity/DrivenCavity.cpp b/apps/gpu/LBM/DrivenCavity/DrivenCavity.cpp
index c02c3b41e633aa6ae63917c6e9e12c2e1a6f2235..42cdce60a11a6d70fc03fe6b734017fefac3eadd 100644
--- a/apps/gpu/LBM/DrivenCavity/DrivenCavity.cpp
+++ b/apps/gpu/LBM/DrivenCavity/DrivenCavity.cpp
@@ -1,375 +1,355 @@
-
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+//      \    \  |    |   ________________________________________________________________
+//       \    \ |    |  |  ______________________________________________________________|
+//        \    \|    |  |  |         __          __     __     __     ______      _______
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file LidDrivenCavity.cpp
+//! \ingroup Applications
+//! \author Martin Schoenherr, Stephan Lenz
+//=======================================================================================
 #define _USE_MATH_DEFINES
-#include <string>
-#include <sstream>
-#include <iostream>
-#include <stdexcept>
-#include <fstream>
 #include <exception>
+#include <fstream>
+#include <iostream>
 #include <memory>
-#include <filesystem>
+#include <sstream>
+#include <stdexcept>
+#include <string>
 
 //////////////////////////////////////////////////////////////////////////
 
 #include "Core/DataTypes.h"
-#include "PointerDefinitions.h"
-
-#include "Core/StringUtilities/StringUtil.h"
-
+#include "Core/LbmOrGks.h"
+#include "Core/Logger/Logger.h"
 #include "Core/VectorTypes.h"
-
-#include "basics/config/ConfigurationFile.h"
-
-#include "logger/Logger.h"
+#include "PointerDefinitions.h"
 
 //////////////////////////////////////////////////////////////////////////
 
+#include "GridGenerator/grid/BoundaryConditions/Side.h"
 #include "GridGenerator/grid/GridBuilder/LevelGridBuilder.h"
 #include "GridGenerator/grid/GridBuilder/MultipleGridBuilder.h"
-#include "GridGenerator/grid/BoundaryConditions/Side.h"
 #include "GridGenerator/grid/GridFactory.h"
 
-#include "GridGenerator/io/SimulationFileWriter/SimulationFileWriter.h"
-#include "GridGenerator/io/GridVTKWriter/GridVTKWriter.h"
-#include "GridGenerator/io/STLReaderWriter/STLReader.h"
-#include "GridGenerator/io/STLReaderWriter/STLWriter.h"
-
 //////////////////////////////////////////////////////////////////////////
 
-#include "VirtualFluids_GPU/LBM/Simulation.h"
+#include "VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h"
 #include "VirtualFluids_GPU/Communication/Communicator.h"
-#include "VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.h"
 #include "VirtualFluids_GPU/DataStructureInitializer/GridProvider.h"
-#include "VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.h"
-#include "VirtualFluids_GPU/Parameter/Parameter.h"
-#include "VirtualFluids_GPU/Output/FileWriter.h"
-#include "VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h"
-
+#include "VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.h"
 #include "VirtualFluids_GPU/GPU/CudaMemoryManager.h"
+#include "VirtualFluids_GPU/LBM/Simulation.h"
+#include "VirtualFluids_GPU/Output/FileWriter.h"
+#include "VirtualFluids_GPU/Parameter/Parameter.h"
 
 //////////////////////////////////////////////////////////////////////////
 
-//#include "GksMeshAdapter/GksMeshAdapter.h"
+// #include "GksMeshAdapter/GksMeshAdapter.h"
+// #include "GksGpu/DataBase/DataBase.h"
+// #include "GksGpu/Initializer/Initializer.h"
+// #include "GksGpu/Parameters/Parameters.h"
+// #include "GksGpu/FlowStateData/FlowStateDataConversion.cuh"
+// #include "GksGpu/BoundaryConditions/BoundaryCondition.h"
+// #include "GksGpu/BoundaryConditions/IsothermalWall.h"
+// #include "GksGpu/TimeStepping/NestedTimeStep.h"
+// #include "GksGpu/Analyzer/ConvergenceAnalyzer.h"
+// #include "GksGpu/Analyzer/CupsAnalyzer.h"
+// #include "GksGpu/CudaUtility/CudaUtility.h"
+// #include "GksGpu/Output/VtkWriter.h"
 
-//#include "GksVtkAdapter/VTKInterface.h"
-//
-//#include "GksGpu/DataBase/DataBase.h"
-//#include "GksGpu/Parameters/Parameters.h"
-//#include "GksGpu/Initializer/Initializer.h"
-//
-//#include "GksGpu/FlowStateData/FlowStateDataConversion.cuh"
-//
-//#include "GksGpu/BoundaryConditions/BoundaryCondition.h"
-//#include "GksGpu/BoundaryConditions/IsothermalWall.h"
-//
-//#include "GksGpu/TimeStepping/NestedTimeStep.h"
-//
-//#include "GksGpu/Analyzer/CupsAnalyzer.h"
-//#include "GksGpu/Analyzer/ConvergenceAnalyzer.h"
-//
-//#include "GksGpu/CudaUtility/CudaUtility.h"
-
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//
-//          U s e r    s e t t i n g s
-//
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-//LbmOrGks lbmOrGks = GKS;
-const LbmOrGks lbmOrGks = LBM;
-
-const real L  = 1.0;
-
-const real Re = 1000.0;
-
-const real velocity  = 1.0;
-
-const real dt = (real)0.5e-3;
-
-const uint nx = 64;
-
-const std::string path("output/");
-const std::string gridPath("grid/");
-
-const std::string simulationName("DrivenCavityChim");
-
-const uint timeStepOut = 10000;
-const uint timeStepEnd = 250000;
-
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
 
-void multipleLevel(const std::string& configPath)
+int main(int argc, char *argv[])
 {
-    logging::Logger::addStream(&std::cout);
-    logging::Logger::setDebugLevel(logging::Logger::Level::INFO_LOW);
-    logging::Logger::timeStamp(logging::Logger::ENABLE);
-    logging::Logger::enablePrintedRankNumbers(logging::Logger::ENABLE);
+    try {
+        //////////////////////////////////////////////////////////////////////////
+        // Simulation parameters
+        //////////////////////////////////////////////////////////////////////////
+        std::string path("./output/DrivenCavity");
+        std::string simulationName("LidDrivenCavity");
 
-    vf::gpu::Communicator& communicator = vf::gpu::Communicator::getInstance();
+        const real L = 1.0;
+        const real Re = 1000.0;
+        const real velocity = 1.0;
+        const real dt = (real)0.5e-3;
+        const uint nx = 64;
 
-    auto gridFactory = GridFactory::make();
-    gridFactory->setTriangularMeshDiscretizationMethod(TriangularMeshDiscretizationMethod::POINT_IN_OBJECT);
-    auto gridBuilder = MultipleGridBuilder::makeShared(gridFactory);
+        const uint timeStepOut = 1000;
+        const uint timeStepEnd = 10000;
 
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        // switch between LBM and GKS solver here
+        // LbmOrGks lbmOrGks = GKS;
+        LbmOrGks lbmOrGks = LBM;
 
-    real dx = L / real(nx);
+        //////////////////////////////////////////////////////////////////////////
+        // setup logger
+        //////////////////////////////////////////////////////////////////////////
 
-    gridBuilder->addCoarseGrid(-0.5 * L, -0.5 * L, -0.5 * L,
-                                0.5 * L,  0.5 * L,  0.5 * L, dx);
+        logging::Logger::addStream(&std::cout);
+        logging::Logger::setDebugLevel(logging::Logger::Level::INFO_LOW);
+        logging::Logger::timeStamp(logging::Logger::ENABLE);
+        logging::Logger::enablePrintedRankNumbers(logging::Logger::ENABLE);
 
-    // gridBuilder->addCoarseGrid(-2.0 * dx, -0.5 * L, -0.5 * L,
-    //                             2.0 * dx,  0.5 * L,  0.5 * L, dx);
+        //////////////////////////////////////////////////////////////////////////
+        // setup gridGenerator
+        //////////////////////////////////////////////////////////////////////////
 
-    auto refBox = new Cuboid(-0.1 * L, -0.1 * L, -0.1 * L,
-                              0.1 * L,  0.1 * L,  0.1 * L);
+        auto gridFactory = GridFactory::make();
+        gridFactory->setTriangularMeshDiscretizationMethod(TriangularMeshDiscretizationMethod::POINT_IN_OBJECT);
+        auto gridBuilder = MultipleGridBuilder::makeShared(gridFactory);
 
-    gridBuilder->addGrid(refBox, 1);
+        //////////////////////////////////////////////////////////////////////////
+        // create grid
+        //////////////////////////////////////////////////////////////////////////
 
-    gridBuilder->setNumberOfLayers(0, 0);
+        real dx = L / real(nx);
 
-    gridBuilder->setPeriodicBoundaryCondition(false, false, false);
+        gridBuilder->addCoarseGrid(-0.5 * L, -0.5 * L, -0.5 * L, 0.5 * L, 0.5 * L, 0.5 * L, dx);
 
-    gridBuilder->buildGrids(lbmOrGks, false); // buildGrids() has to be called before setting the BCs!!!!
+        gridBuilder->setPeriodicBoundaryCondition(false, false, false);
 
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        gridBuilder->buildGrids(lbmOrGks, false);
 
-    if( lbmOrGks == LBM )
-    {
+        //////////////////////////////////////////////////////////////////////////
+        // branch between LBM and GKS
+        //////////////////////////////////////////////////////////////////////////
 
-        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-        vf::basics::ConfigurationFile config;
-        config.load(configPath);
+        if (lbmOrGks == LBM) {
+            SPtr<Parameter> para = std::make_shared<Parameter>();
+            BoundaryConditionFactory bcFactory = BoundaryConditionFactory();
+            vf::gpu::Communicator &communicator = vf::gpu::Communicator::getInstance();
 
-        SPtr<Parameter> para = std::make_shared<Parameter>(config, communicator.getNummberOfProcess(), communicator.getPID());
-        BoundaryConditionFactory bcFactory = BoundaryConditionFactory();
+            //////////////////////////////////////////////////////////////////////////
+            // compute parameters in lattice units
+            //////////////////////////////////////////////////////////////////////////
 
-        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+            const real velocityLB = velocity * dt / dx; // LB units
 
-        const real velocityLB = velocity * dt / dx; // LB units
+            const real vxLB = velocityLB / sqrt(2.0); // LB units
+            const real vyLB = velocityLB / sqrt(2.0); // LB units
 
-        const real vx = velocityLB / (real)sqrt(2.0); // LB units
-        const real vy = velocityLB / (real)sqrt(2.0); // LB units
+            const real viscosityLB = nx * velocityLB / Re; // LB units
 
-        const real viscosityLB = nx * velocityLB / Re; // LB units
+            *logging::out << logging::Logger::INFO_HIGH << "velocity  [dx/dt] = " << velocityLB << " \n";
+            *logging::out << logging::Logger::INFO_HIGH << "viscosity [dx^2/dt] = " << viscosityLB << "\n";
 
-        VF_LOG_INFO("velocity  [dx/dt] = {}", velocityLB);
-        VF_LOG_INFO("viscosity [dx^2/dt] = {}", viscosityLB);
+            //////////////////////////////////////////////////////////////////////////
+            // set parameters
+            //////////////////////////////////////////////////////////////////////////
 
-        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+            para->setOutputPath(path);
+            para->setOutputPrefix(simulationName);
 
-        para->setDevices(std::vector<uint>{(uint)0});
+            para->setPrintFiles(true);
 
-        para->setOutputPath( path ); // optional, default is output/
-        para ->setGridPath( gridPath );  // optional, default is grid/
+            para->setVelocityLB(velocityLB);
+            para->setViscosityLB(viscosityLB);
 
-        para->setOutputPrefix( simulationName );
+            para->setVelocityRatio(velocity / velocityLB);
+            para->setDensityRatio(1.0);
 
-        para->setPrintFiles(true);
+            para->setTimestepOut(timeStepOut);
+            para->setTimestepEnd(timeStepEnd);
 
-        para->setMaxLevel(2);
+            para->setMainKernel("CumulantK17CompChim");
 
-        para->setVelocity(velocityLB);
-        para->setViscosity(viscosityLB);
+            //////////////////////////////////////////////////////////////////////////
+            // set boundary conditions
+            //////////////////////////////////////////////////////////////////////////
 
-        para->setVelocityRatio(velocity / velocityLB);
+            gridBuilder->setNoSlipBoundaryCondition(SideType::PX);
+            gridBuilder->setNoSlipBoundaryCondition(SideType::MX);
+            gridBuilder->setNoSlipBoundaryCondition(SideType::PY);
+            gridBuilder->setNoSlipBoundaryCondition(SideType::MY);
+            gridBuilder->setVelocityBoundaryCondition(SideType::PZ, vxLB, vyLB, 0.0);
+            gridBuilder->setNoSlipBoundaryCondition(SideType::MZ);
 
-        //para->setMainKernel("CumulantK17CompChim");
+            bcFactory.setNoSlipBoundaryCondition(BoundaryConditionFactory::NoSlipBC::NoSlipBounceBack);
+            bcFactory.setVelocityBoundaryCondition(BoundaryConditionFactory::VelocityBC::VelocitySimpleBounceBackCompressible);
 
-        para->setInitialCondition([&](real coordX, real coordY, real coordZ, real &rho, real &vx, real &vy, real &vz) {
-            rho = (real)0.0;
-            vx  = (real)0.0; //(6 * velocityLB * coordZ * (L - coordZ) / (L * L));
-            vy  = (real)0.0;
-            vz  = (real)0.0;
-        });
+            //////////////////////////////////////////////////////////////////////////
+            // set copy mesh to simulation
+            //////////////////////////////////////////////////////////////////////////
 
-        para->setTOut( timeStepOut );
-        para->setTEnd( timeStepEnd );
+            auto cudaMemoryManager = std::make_shared<CudaMemoryManager>(para);
+            SPtr<GridProvider> gridGenerator =
+                GridProvider::makeGridGenerator(gridBuilder, para, cudaMemoryManager, communicator);
 
-        /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+            para->setCalcParticles(false);
+            para->setCalcMedian(false);
+            para->setCalcTurbulenceIntensity(false);
+            para->setDoCheckPoint(false);
+            para->setUseMeasurePoints(false);
+            para->setDiffOn(false);
+            para->setCalcPlaneConc(false);
+            para->setUseWale(false);
+            para->setSimulatePorousMedia(false);
+            //////////////////////////////////////////////////////////////////////////
+            // run simulation
+            //////////////////////////////////////////////////////////////////////////
 
-        gridBuilder->setNoSlipBoundaryCondition(SideType::PX);
-        gridBuilder->setNoSlipBoundaryCondition(SideType::MX);
-        gridBuilder->setNoSlipBoundaryCondition(SideType::PY);
-        gridBuilder->setNoSlipBoundaryCondition(SideType::MY);
-        gridBuilder->setVelocityBoundaryCondition(SideType::PZ, vx, vx, 0.0);
-        gridBuilder->setNoSlipBoundaryCondition(SideType::MZ);
+            Simulation sim(para, cudaMemoryManager, communicator, *gridGenerator, &bcFactory);
+            sim.run();
+        } // else {
+        //     CudaUtility::setCudaDevice(0);
 
-        bcFactory.setNoSlipBoundaryCondition(BoundaryConditionFactory::NoSlipBC::NoSlipBounceBack);
-        bcFactory.setVelocityBoundaryCondition(BoundaryConditionFactory::VelocityBC::VelocitySimpleBounceBackCompressible);
+        //     Parameters parameters;
 
-        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //     //////////////////////////////////////////////////////////////////////////
+        //     // compute remaining parameters
+        //     //////////////////////////////////////////////////////////////////////////
 
-        gridBuilder->writeGridsToVtk(para->getGridPath());
+        //     const real vx = velocity / sqrt(2.0);
+        //     const real vy = velocity / sqrt(2.0);
 
-        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //     parameters.K  = 2.0;
+        //     parameters.Pr = 1.0;
 
-        auto cudaMemoryManager = std::make_shared<CudaMemoryManager>(para);
+        //     const real Ma = (real)0.1;
 
-        auto gridGenerator = GridProvider::makeGridGenerator(gridBuilder, para, cudaMemoryManager, communicator);
+        //     real rho = 1.0;
 
-        Simulation sim(para, cudaMemoryManager, communicator, *gridGenerator, &bcFactory);
-        sim.run();
-
-        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    }
-    else
-    {
-     //   CudaUtility::setCudaDevice(0);
-     //
-     //   Parameters parameters;
+        //     real cs     = velocity / Ma;
+        //     real lambda = c1o2 * ((parameters.K + 5.0) / (parameters.K + 3.0)) / (cs * cs);
 
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //     const real mu = velocity * L * rho / Re;
 
-        //const real vx = velocity / sqrt(2.0);
-        //const real vy = velocity / sqrt(2.0);
+        //     *logging::out << logging::Logger::INFO_HIGH << "mu  = " << mu << " m^2/s\n";
 
-     //   parameters.K  = 2.0;
-     //   parameters.Pr = 1.0;
-     //
-     //   const real Ma = 0.1;
+        //     *logging::out << logging::Logger::INFO_HIGH << "CFL = " << dt * (velocity + cs) / dx << "\n";
 
-     //   real rho = 1.0;
+        //     //////////////////////////////////////////////////////////////////////////
+        //     // set parameters
+        //     //////////////////////////////////////////////////////////////////////////
 
-     //   real cs = velocity / Ma;
-     //   real lambda = c1o2 * ( ( parameters.K + 5.0 ) / ( parameters.K + 3.0 ) ) / ( cs * cs );
+        //     parameters.mu = mu;
 
-     //   const real mu = velocity * L * rho / Re;
+        //     parameters.dt = dt;
+        //     parameters.dx = dx;
 
-     //   *logging::out << logging::Logger::INFO_HIGH << "mu  = " << mu << " m^2/s\n";
+        //     parameters.lambdaRef = lambda;
 
-     //   *logging::out << logging::Logger::INFO_HIGH << "CFL = " << dt * ( velocity + cs ) / dx << "\n";
+        //     //////////////////////////////////////////////////////////////////////////
+        //     // set copy mesh to simulation
+        //     //////////////////////////////////////////////////////////////////////////
 
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //     GksMeshAdapter meshAdapter(gridBuilder);
 
-     //   parameters.mu = mu;
+        //     meshAdapter.inputGrid();
 
-     //   parameters.dt = dt;
-     //   parameters.dx = dx;
+        //     auto dataBase = std::make_shared<DataBase>("GPU");
 
-     //   parameters.lambdaRef = lambda;
+        //     //////////////////////////////////////////////////////////////////////////
+        //     // set boundary conditions
+        //     //////////////////////////////////////////////////////////////////////////
 
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //     SPtr<BoundaryCondition> bcLid =
+        //         std::make_shared<IsothermalWall>(dataBase, Vec3(vx, vy, 0.0), lambda, false);
+        //     SPtr<BoundaryCondition> bcWall =
+        //         std::make_shared<IsothermalWall>(dataBase, Vec3(0.0, 0.0, 0.0), lambda, false);
 
-     //   GksMeshAdapter meshAdapter( gridBuilder );
+        //     bcLid->findBoundaryCells(meshAdapter, false, [&](Vec3 center) {
+        //         return center.z > 0.5 && center.x > -0.5 && center.x < 0.5 && center.y > -0.5 && center.y < 0.5;
+        //     });
 
-     //   meshAdapter.inputGrid();
+        //     bcWall->findBoundaryCells(meshAdapter, true, [&](Vec3 center) {
+        //         return center.x < -0.5 || center.x > 0.5 || center.y < -0.5 || center.y > 0.5 || center.z < -0.5;
+        //     });
 
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //     dataBase->boundaryConditions.push_back(bcLid);
+        //     dataBase->boundaryConditions.push_back(bcWall);
 
-     //   auto dataBase = std::make_shared<DataBase>( "GPU" );
+        //     //////////////////////////////////////////////////////////////////////////
+        //     // set initial condition and upload mesh and initial condition to GPGPU
+        //     //////////////////////////////////////////////////////////////////////////
 
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //     dataBase->setMesh(meshAdapter);
 
-     //   SPtr<BoundaryCondition> bcLid  = std::make_shared<IsothermalWall>( dataBase, Vec3(  vx,  vy, 0.0 ), lambda, false );
-     //   SPtr<BoundaryCondition> bcWall = std::make_shared<IsothermalWall>( dataBase, Vec3( 0.0, 0.0, 0.0 ), lambda, false );
+        //     Initializer::interpret(dataBase, [&](Vec3 cellCenter) -> ConservedVariables {
+        //         return toConservedVariables(PrimitiveVariables(rho, 0.0, 0.0, 0.0, lambda), parameters.K);
+        //     });
 
-     //   bcLid->findBoundaryCells ( meshAdapter, true,  [&](Vec3 center){ return center.z > 0.5; } );
-     //   bcWall->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.z < 0.5; } );
+        //     dataBase->copyDataHostToDevice();
 
-     //   dataBase->boundaryConditions.push_back( bcLid  );
-     //   dataBase->boundaryConditions.push_back( bcWall );
+        //     Initializer::initializeDataUpdate(dataBase);
 
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //     VtkWriter::write(dataBase, parameters, path + "/" + simulationName + "_0");
 
-     //   dataBase->setMesh( meshAdapter );
+        //     //////////////////////////////////////////////////////////////////////////
+        //     // set analyzers
+        //     //////////////////////////////////////////////////////////////////////////
 
-     //   Initializer::interpret(dataBase, [&] ( Vec3 cellCenter ) -> ConservedVariables {
+        //     CupsAnalyzer cupsAnalyzer(dataBase, false, 60.0, true, 10000);
 
-     //       return toConservedVariables( PrimitiveVariables( rho, 0.0, 0.0, 0.0, lambda ), parameters.K );
-     //   });
+        //     ConvergenceAnalyzer convergenceAnalyzer(dataBase, 10000);
 
-     //   dataBase->copyDataHostToDevice();
+        //     cupsAnalyzer.start();
 
-     //   Initializer::initializeDataUpdate(dataBase);
+        //     //////////////////////////////////////////////////////////////////////////
+        //     // run simulation
+        //     //////////////////////////////////////////////////////////////////////////
 
-     //   writeVtkXML( dataBase, parameters, 0, path + simulationName + "_0" );
+        //     for (uint iter = 1; iter <= timeStepEnd; iter++) {
+        //         TimeStepping::nestedTimeStep(dataBase, parameters, 0);
 
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //         if (iter % timeStepOut == 0) {
+        //             dataBase->copyDataDeviceToHost();
 
-     //   CupsAnalyzer cupsAnalyzer( dataBase, false, 60.0, true, 10000 );
+        //             VtkWriter::write(dataBase, parameters, path + "/" + simulationName + "_" + std::to_string(iter));
+        //         }
 
-     //   ConvergenceAnalyzer convergenceAnalyzer( dataBase, 10000 );
+        //         int crashCellIndex = dataBase->getCrashCellIndex();
+        //         if (crashCellIndex >= 0) {
+        //             *logging::out << logging::Logger::LOGGER_ERROR
+        //                           << "Simulation crashed at CellIndex = " << crashCellIndex << "\n";
+        //             dataBase->copyDataDeviceToHost();
+        //             VtkWriter::write(dataBase, parameters, path + "/" + simulationName + "_" + std::to_string(iter));
 
-     //   cupsAnalyzer.start();
+        //             break;
+        //         }
 
-     //   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //         dataBase->getCrashCellIndex();
 
-     //   for( uint iter = 1; iter <= timeStepEnd; iter++ )
-     //   {
-     //       TimeStepping::nestedTimeStep(dataBase, parameters, 0);
+        //         cupsAnalyzer.run(iter, parameters.dt);
 
-     //       if( iter % timeStepOut == 0 )
-     //       {
-     //           dataBase->copyDataDeviceToHost();
+        //         convergenceAnalyzer.run(iter);
+        //     }
+        // }
+    } catch (const std::bad_alloc &e) {
 
-     //           writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( iter ) );
-     //       }
-     //
-     //       int crashCellIndex = dataBase->getCrashCellIndex();
-     //       if( crashCellIndex >= 0 )
-     //       {
-     //           *logging::out << logging::Logger::LOGGER_ERROR << "Simulation Crashed at CellIndex = " << crashCellIndex << "\n";
-     //           dataBase->copyDataDeviceToHost();
-     //           writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( iter ) );
+        *logging::out << logging::Logger::LOGGER_ERROR << "Bad Alloc:" << e.what() << "\n";
+    } catch (const std::exception &e) {
 
-     //           break;
-     //       }
+        *logging::out << logging::Logger::LOGGER_ERROR << e.what() << "\n";
+    } catch (std::string &s) {
 
-     //       dataBase->getCrashCellIndex();
-
-     //       cupsAnalyzer.run( iter, parameters.dt );
-
-     //       convergenceAnalyzer.run( iter );
-     //   }
-    }
-}
-
-int main( int argc, char* argv[])
-{
-    try
-    {
-        vf::logging::Logger::initalizeLogger();
-
-        // assuming that the config files is stored parallel to this file.
-        std::filesystem::path filePath = __FILE__;
-        filePath.replace_filename("configDrivenCavity.txt");
-
-        multipleLevel(filePath.string());
-    }
-    catch (const spdlog::spdlog_ex &ex) {
-        std::cout << "Log initialization failed: " << ex.what() << std::endl;
-    }
-    catch (const std::bad_alloc& e)
-    {
-        VF_LOG_CRITICAL("Bad Alloc: {}", e.what());
-    }
-    catch (const std::exception& e)
-    {
-        VF_LOG_CRITICAL("exception: {}", e.what());
-    }
-    catch (...)
-    {
-        VF_LOG_CRITICAL("Unknown exception!");
+        *logging::out << logging::Logger::LOGGER_ERROR << s << "\n";
+    } catch (...) {
+        *logging::out << logging::Logger::LOGGER_ERROR << "Unknown exception!\n";
     }
 
-   return 0;
+    return 0;
 }
diff --git a/apps/gpu/LBM/DrivenCavityMultiGPU/DrivenCavityMultiGPU.cpp b/apps/gpu/LBM/DrivenCavityMultiGPU/DrivenCavityMultiGPU.cpp
index 0376e376b365c8c8282cc3cd79a5c365b416193f..c62f42f9943ad17a4b5f6c61838c2b07ab5f7c15 100644
--- a/apps/gpu/LBM/DrivenCavityMultiGPU/DrivenCavityMultiGPU.cpp
+++ b/apps/gpu/LBM/DrivenCavityMultiGPU/DrivenCavityMultiGPU.cpp
@@ -137,8 +137,8 @@ void multipleLevel(const std::string &configPath)
     *logging::out << logging::Logger::INFO_HIGH << "velocity  [dx/dt] = " << velocityLB << " \n";
     *logging::out << logging::Logger::INFO_HIGH << "viscosity [dx^2/dt] = " << viscosityLB << "\n";
 
-    para->setVelocity(velocityLB);
-    para->setViscosity(viscosityLB);
+    para->setVelocityLB(velocityLB);
+    para->setViscosityLB(viscosityLB);
     para->setVelocityRatio(velocity / velocityLB);
     para->setDensityRatio((real)1.0); // correct value?
 
diff --git a/apps/gpu/LBM/MusselOyster/MusselOyster.cpp b/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
index 0bef3e75ff93e42f84477431bf78e7cbf7702b17..7ddffb06272d6c650762b87afc181af6ecca51da 100644
--- a/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
+++ b/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
@@ -140,8 +140,8 @@ void multipleLevel(const std::string &configPath)
     real referenceLength = 1.0 / dxGrid; // heightBivalve / dxGrid
     real viscosityLB     = (vxLB * referenceLength) / Re;
 
-    para->setVelocity(vxLB);
-    para->setViscosity(viscosityLB);
+    para->setVelocityLB(vxLB);
+    para->setViscosityLB(viscosityLB);
     para->setVelocityRatio((real)58.82352941);
     para->setViscosityRatio((real)0.058823529);
     para->setDensityRatio((real)998.0);
diff --git a/apps/gpu/LBM/SphereGPU/Sphere.cpp b/apps/gpu/LBM/SphereGPU/Sphere.cpp
index f0d629afecca44f917d002543af7689e0889f42b..a4a6a56ba594b1b933fa971356f57db91fb2adf6 100644
--- a/apps/gpu/LBM/SphereGPU/Sphere.cpp
+++ b/apps/gpu/LBM/SphereGPU/Sphere.cpp
@@ -1,4 +1,35 @@
-
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+//      \    \  |    |   ________________________________________________________________
+//       \    \ |    |  |  ______________________________________________________________|
+//        \    \|    |  |  |         __          __     __     __     ______      _______
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file LidDrivenCavity.cpp
+//! \ingroup Applications
+//! \author Martin Schoenherr, Stephan Lenz, Anna Wellmann
+//=======================================================================================
 #define _USE_MATH_DEFINES
 #include <exception>
 #include <filesystem>
@@ -13,15 +44,11 @@
 //////////////////////////////////////////////////////////////////////////
 
 #include "Core/DataTypes.h"
-#include "PointerDefinitions.h"
-
-#include "Core/StringUtilities/StringUtil.h"
-
+#include "Core/LbmOrGks.h"
+#include "Core/Logger/Logger.h"
 #include "Core/VectorTypes.h"
-
-#include <basics/config/ConfigurationFile.h>
-
-#include <logger/Logger.h>
+#include "PointerDefinitions.h"
+#include "config/ConfigurationFile.h"
 
 //////////////////////////////////////////////////////////////////////////
 
@@ -30,187 +57,179 @@
 #include "GridGenerator/grid/GridBuilder/MultipleGridBuilder.h"
 #include "GridGenerator/grid/GridFactory.h"
 
-#include "GridGenerator/io/GridVTKWriter/GridVTKWriter.h"
-#include "GridGenerator/io/STLReaderWriter/STLReader.h"
-#include "GridGenerator/io/STLReaderWriter/STLWriter.h"
-#include "GridGenerator/io/SimulationFileWriter/SimulationFileWriter.h"
-
 #include "GridGenerator/geometries/Sphere/Sphere.h"
+#include "GridGenerator/geometries/TriangularMesh/TriangularMesh.h"
 
 //////////////////////////////////////////////////////////////////////////
 
-#include "VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h"
-#include "VirtualFluids_GPU/Communication/Communicator.h"
 #include "VirtualFluids_GPU/DataStructureInitializer/GridProvider.h"
-#include "VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.h"
 #include "VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.h"
+#include "VirtualFluids_GPU/GPU/CudaMemoryManager.h"
+#include "VirtualFluids_GPU/Communication/Communicator.h"
 #include "VirtualFluids_GPU/LBM/Simulation.h"
 #include "VirtualFluids_GPU/Output/FileWriter.h"
 #include "VirtualFluids_GPU/Parameter/Parameter.h"
+#include "VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h"
+#include "VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h"
 
-#include "VirtualFluids_GPU/GPU/CudaMemoryManager.h"
-
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//
-//          U s e r    s e t t i n g s
-//
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-const real L = 1.0;
-const real dSphere = 0.2;
-const real Re = 1000.0; // related to the sphere's diameter
-const real velocity = 1.0;
-const real dt = (real)0.5e-3;
-const uint nx = 64;
-
-const uint timeStepOut = 10000;
-const uint timeStepEnd = 100000;
-
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
 
-void multipleLevel(const std::string &configPath)
+int main(int argc, char *argv[])
 {
-    logging::Logger::addStream(&std::cout);
-    logging::Logger::setDebugLevel(logging::Logger::Level::INFO_LOW);
-    logging::Logger::timeStamp(logging::Logger::ENABLE);
-    logging::Logger::enablePrintedRankNumbers(logging::Logger::ENABLE);
-
-    vf::gpu::Communicator &communicator = vf::gpu::Communicator::getInstance();
-
-    auto gridFactory = GridFactory::make();
-    gridFactory->setTriangularMeshDiscretizationMethod(TriangularMeshDiscretizationMethod::POINT_IN_OBJECT);
-    auto gridBuilder = MultipleGridBuilder::makeShared(gridFactory);
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    real dx = L / real(nx);
-
-    gridBuilder->addCoarseGrid(-1.0 * L, -0.5 * L, -0.5 * L, 
-                                6.0 * L,  0.5 * L,  0.5 * L, dx);
-
-    Object *sphere = new Sphere(0.0, 0.0, 0.0, dSphere / 2.0);
-    gridBuilder->addGeometry(sphere);
-
-    gridBuilder->setPeriodicBoundaryCondition(false, false, false);
-
-    gridBuilder->buildGrids(LBM, false); // buildGrids() has to be called before setting the BCs!!!!
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    vf::basics::ConfigurationFile config;
-    config.load(configPath);
-
-    SPtr<Parameter> para =
-        std::make_shared<Parameter>(config, communicator.getNummberOfProcess(), communicator.getPID());
-    BoundaryConditionFactory bcFactory = BoundaryConditionFactory();
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    const real velocityLB = velocity * dt / dx; // LB units
-
-    const real viscosityLB = (dSphere / dx) * velocityLB / Re; // LB units
-
-    VF_LOG_INFO("velocity  [dx/dt] = {}", velocityLB);
-    VF_LOG_INFO("viscosity [dx^2/dt] = {}", viscosityLB);
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    para->setDevices(std::vector<uint>{ (uint)0 });
-
-    para->setFName(para->getOutputPath() + "/" + para->getOutputPrefix());
-
-    para->setPrintFiles(true);
-
-    para->setMaxLevel(2);
-
-    para->setVelocity(velocityLB);
-    para->setViscosity(viscosityLB);
-
-    para->setVelocityRatio(velocity / velocityLB);
-    para->setDensityRatio((real)1.0);
-
-    // para->setMainKernel("CumulantK17CompChim");
-
-    // para->setInitialCondition([&](real coordX, real coordY, real coordZ, real &rho, real &vx, real &vy, real &vz) {
-    //     rho = (real)0.0;
-    //     vx =  (real)velocityLB;
-    //     vy =  (real)0.0;
-    //     vz =  (real)0.0;
-    // });
-
-    para->setTOut(timeStepOut);
-    para->setTEnd(timeStepEnd);
-
-    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    gridBuilder->setVelocityBoundaryCondition(SideType::MX, velocityLB, 0.0, 0.0);
-
-    gridBuilder->setSlipBoundaryCondition(SideType::PY, 0.0, 0.0, 0.0);
-    gridBuilder->setSlipBoundaryCondition(SideType::MY, 0.0, 0.0, 0.0);
-    gridBuilder->setSlipBoundaryCondition(SideType::PZ, 0.0, 0.0, 0.0);
-    gridBuilder->setSlipBoundaryCondition(SideType::MZ, 0.0, 0.0, 0.0);
-
-    // gridBuilder->setNoSlipBoundaryCondition(SideType::GEOMETRY); // not working yet, use veloBC
-    gridBuilder->setVelocityBoundaryCondition(SideType::GEOMETRY, 0.0, velocityLB * 2.0, 0.0);
-    gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0); // set pressure boundary condition last
-
-    bcFactory.setNoSlipBoundaryCondition(BoundaryConditionFactory::NoSlipBC::NoSlipCompressible);
-    bcFactory.setVelocityBoundaryCondition(BoundaryConditionFactory::VelocityBC::VelocityCompressible);
-    bcFactory.setSlipBoundaryCondition(BoundaryConditionFactory::SlipBC::SlipCompressible);
-    bcFactory.setPressureBoundaryCondition(BoundaryConditionFactory::PressureBC::PressureNonEquilibriumCompressible);
-    bcFactory.setGeometryBoundaryCondition(BoundaryConditionFactory::NoSlipBC::NoSlipCompressible);
-    // bcFactory.setGeometryBoundaryCondition(BoundaryConditionFactory::VelocityBC::VelocityCompressible);
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    // gridBuilder->writeGridsToVtk("grid/");
+    try {
+        //////////////////////////////////////////////////////////////////////////
+        // Simulation parameters
+        //////////////////////////////////////////////////////////////////////////
+
+        const bool useConfigFile = true;
+
+        const real L = 1.0;
+        const real dSphere = 0.2;
+        const real Re = 1000.0; // related to the sphere's diameter
+        const real velocity = 1.0;
+        const real dt = (real)0.5e-3;
+        const uint nx = 64;
+
+        const uint timeStepOut = 10000;
+        const uint timeStepEnd = 100000;
+
+        //////////////////////////////////////////////////////////////////////////
+        // setup logger
+        //////////////////////////////////////////////////////////////////////////
+
+        logging::Logger::addStream(&std::cout);
+        logging::Logger::setDebugLevel(logging::Logger::Level::INFO_LOW);
+        logging::Logger::timeStamp(logging::Logger::ENABLE);
+        logging::Logger::enablePrintedRankNumbers(logging::Logger::ENABLE);
+
+        //////////////////////////////////////////////////////////////////////////
+        // setup simulation parameters (with or without config file)
+        //////////////////////////////////////////////////////////////////////////
+
+        vf::gpu::Communicator& communicator = vf::gpu::Communicator::getInstance();;
+        SPtr<Parameter> para;
+        BoundaryConditionFactory bcFactory = BoundaryConditionFactory();
+        vf::basics::ConfigurationFile config;
+        if (useConfigFile) {
+            //////////////////////////////////////////////////////////////////////////
+            // read simulation parameters from config file
+            //////////////////////////////////////////////////////////////////////////
+
+            // assuming that a config files is stored parallel to this file.
+            std::filesystem::path configPath = __FILE__;
+
+            // the config file's default name can be replaced by passing a command line argument
+            std::string configName("config.txt");
+            if (argc == 2) {
+                configName = argv[1];
+                std::cout << "Using configFile command line argument: " << configName << std::endl;
+            }
+
+            configPath.replace_filename(configName);
+            config.load(configPath.string());
+
+            para = std::make_shared<Parameter>(config);
+        } else {
+            para = std::make_shared<Parameter>();
+        }
+
+        //////////////////////////////////////////////////////////////////////////
+        // setup gridGenerator
+        //////////////////////////////////////////////////////////////////////////
+
+        auto gridFactory = GridFactory::make();
+        gridFactory->setTriangularMeshDiscretizationMethod(TriangularMeshDiscretizationMethod::POINT_IN_OBJECT);
+        auto gridBuilder = MultipleGridBuilder::makeShared(gridFactory);
+
+        //////////////////////////////////////////////////////////////////////////
+        // create grid
+        //////////////////////////////////////////////////////////////////////////
+
+        real dx = L / real(nx);
+        gridBuilder->addCoarseGrid(-1.0 * L, -0.8 * L, -0.8 * L,
+                                    6.0 * L,  0.8 * L,  0.8 * L, dx);
+
+        // use primitive
+        Object *sphere = new Sphere(0.0, 0.0, 0.0, dSphere / 2.0);
+
+        // use stl
+        // std::string stlPath = "stl/sphere02.stl";
+        // if (useConfigFile && config.contains("STLPath")) {
+        //     stlPath = config.getValue<std::string>("STLPath");
+        // }
+        // std::cout << "Reading stl from " << stlPath << "." << std::endl;
+        // Object *sphere = TriangularMesh::make(stlPath);
+
+        gridBuilder->addGeometry(sphere);
+        gridBuilder->setPeriodicBoundaryCondition(false, false, false);
+        gridBuilder->buildGrids(LBM, false);  // buildGrids() has to be called before setting the BCs!!!!
+
+        //////////////////////////////////////////////////////////////////////////
+        // compute parameters in lattice units
+        //////////////////////////////////////////////////////////////////////////
+
+        const real velocityLB = velocity * dt / dx; // LB units
+        const real viscosityLB =  (dSphere / dx) * velocityLB / Re; // LB units
+
+        *logging::out << logging::Logger::INFO_HIGH << "velocity  [dx/dt] = " << velocityLB << " \n";
+        *logging::out << logging::Logger::INFO_HIGH << "viscosity [dx^2/dt] = " << viscosityLB << "\n";
+
+        //////////////////////////////////////////////////////////////////////////
+        // set parameters
+        //////////////////////////////////////////////////////////////////////////
+
+        para->setPrintFiles(true);
 
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        para->setVelocityLB(velocityLB);
+        para->setViscosityLB(viscosityLB);
+
+        para->setVelocityRatio(velocity / velocityLB);
+        para->setDensityRatio((real)1.0);
+
+        para->setTimestepOut(timeStepOut);
+        para->setTimestepEnd(timeStepEnd);
 
-    auto cudaMemoryManager = std::make_shared<CudaMemoryManager>(para);
+        //////////////////////////////////////////////////////////////////////////
+        // set boundary conditions
+        //////////////////////////////////////////////////////////////////////////
 
-    auto gridGenerator = GridProvider::makeGridGenerator(gridBuilder, para, cudaMemoryManager, communicator);
+        gridBuilder->setVelocityBoundaryCondition(SideType::MX, velocityLB, 0.0, 0.0);
+
+        gridBuilder->setSlipBoundaryCondition(SideType::PY, 0.0, 0.0, 0.0);
+        gridBuilder->setSlipBoundaryCondition(SideType::MY, 0.0, 0.0, 0.0);
+        gridBuilder->setSlipBoundaryCondition(SideType::PZ, 0.0, 0.0, 0.0);
+        gridBuilder->setSlipBoundaryCondition(SideType::MZ, 0.0, 0.0, 0.0);
 
-    Simulation sim(para, cudaMemoryManager, communicator, *gridGenerator, &bcFactory);
-    sim.run();
+        gridBuilder->setVelocityBoundaryCondition(SideType::GEOMETRY, 0.0, 0.0, 0.0);
+        gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0); // set pressure boundary condition last
+
+        bcFactory.setVelocityBoundaryCondition(BoundaryConditionFactory::VelocityBC::VelocityCompressible);
+        bcFactory.setSlipBoundaryCondition(BoundaryConditionFactory::SlipBC::SlipCompressible);
+        bcFactory.setPressureBoundaryCondition(BoundaryConditionFactory::PressureBC::PressureNonEquilibriumCompressible);
+        bcFactory.setGeometryBoundaryCondition(BoundaryConditionFactory::NoSlipBC::NoSlipCompressible);
+        
+        //////////////////////////////////////////////////////////////////////////
+        // setup to copy mesh to simulation
+        //////////////////////////////////////////////////////////////////////////
 
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-}
+        auto cudaMemoryManager = std::make_shared<CudaMemoryManager>(para);
+        SPtr<GridProvider> gridGenerator = GridProvider::makeGridGenerator(gridBuilder, para, cudaMemoryManager, communicator);
 
-int main(int argc, char *argv[])
-{
-    try {
-        vf::logging::Logger::initalizeLogger();
+        //////////////////////////////////////////////////////////////////////////
+        // run simulation
+        //////////////////////////////////////////////////////////////////////////
 
-        // assuming that the config files is stored parallel to this file.
-        std::filesystem::path filePath = __FILE__;
-        filePath.replace_filename("config.txt");
+        Simulation sim(para, cudaMemoryManager, communicator, *gridGenerator, &bcFactory);
+        sim.run();
 
-        multipleLevel(filePath.string());
-    } catch (const spdlog::spdlog_ex &ex) {
-        std::cout << "Log initialization failed: " << ex.what() << std::endl;
     } catch (const std::bad_alloc &e) {
-        VF_LOG_CRITICAL("Bad Alloc: {}", e.what());
+        *logging::out << logging::Logger::LOGGER_ERROR << "Bad Alloc:" << e.what() << "\n";
     } catch (const std::exception &e) {
-        VF_LOG_CRITICAL("exception: {}", e.what());
+        *logging::out << logging::Logger::LOGGER_ERROR << e.what() << "\n";
+    } catch (std::string &s) {
+        *logging::out << logging::Logger::LOGGER_ERROR << s << "\n";
     } catch (...) {
-        VF_LOG_CRITICAL("Unknown exception!");
+        *logging::out << logging::Logger::LOGGER_ERROR << "Unknown exception!\n";
     }
 
     return 0;
diff --git a/apps/gpu/LBM/SphereScaling/SphereScaling.cpp b/apps/gpu/LBM/SphereScaling/SphereScaling.cpp
index e9058e455e9d41173e3c4d500a456cb4f71d6502..250414bbf0efbe09d24471309f7d6aaa280c5477 100644
--- a/apps/gpu/LBM/SphereScaling/SphereScaling.cpp
+++ b/apps/gpu/LBM/SphereScaling/SphereScaling.cpp
@@ -133,8 +133,8 @@ void multipleLevel(const std::string &configPath)
     real vxLB        = (real)0.0005; // LB units
     real viscosityLB = 0.001;        //(vxLB * dxGrid) / Re;
 
-    para->setVelocity(vxLB);
-    para->setViscosity(viscosityLB);
+    para->setVelocityLB(vxLB);
+    para->setViscosityLB(viscosityLB);
     para->setVelocityRatio((real)58.82352941);
     para->setViscosityRatio((real)0.058823529);
     para->setDensityRatio((real)998.0);
diff --git a/apps/gpu/LBM/TGV_3D/TGV_3D.cpp b/apps/gpu/LBM/TGV_3D/TGV_3D.cpp
index 34793f04d901ac9a412c8085f169e9e79c61c3e4..904f1e29d11a940b41ec6b413ea3545638089536 100644
--- a/apps/gpu/LBM/TGV_3D/TGV_3D.cpp
+++ b/apps/gpu/LBM/TGV_3D/TGV_3D.cpp
@@ -217,9 +217,9 @@ void multipleLevel(const std::string& configPath)
     para->setTEnd( 40 * lround(L/velocity) );	
 	para->setTOut(  5 * lround(L/velocity) );
 
-    para->setVelocity( velocity );
+    para->setVelocityLB( velocity );
 
-    para->setViscosity( viscosity );
+    para->setViscosityLB( viscosity );
 
     para->setVelocityRatio( 1.0 / velocity );
 
diff --git a/apps/gpu/LBM/TGV_3D_MultiGPU/TGV_3D_MultiGPU.cpp b/apps/gpu/LBM/TGV_3D_MultiGPU/TGV_3D_MultiGPU.cpp
index 3c2362c1c997ab7eb226d7aa9c51b40b1e7e5605..9488e1e7671eac1e48dfe9e1ef85db878df65cf5 100644
--- a/apps/gpu/LBM/TGV_3D_MultiGPU/TGV_3D_MultiGPU.cpp
+++ b/apps/gpu/LBM/TGV_3D_MultiGPU/TGV_3D_MultiGPU.cpp
@@ -268,9 +268,9 @@ void multipleLevel(const std::string& configPath)
     para->setTEnd( 1000 );	
 	//para->setTOut(    1 );
 
-    para->setVelocity( velocity );
+    para->setVelocityLB( velocity );
 
-    para->setViscosity( viscosity );
+    para->setViscosityLB( viscosity );
 
     para->setVelocityRatio( 1.0 / velocity );
 
diff --git a/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp b/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp
index 72577799763e607d1cc5bea0d48e807594b9a5f8..b06cda9399fbbca78131a8cf0fb8753b1556aff0 100644
--- a/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp
+++ b/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp
@@ -225,8 +225,8 @@ void multipleLevel(const std::string& configPath)
 
     para->setMaxLevel(maxLevel);
 
-    para->setVelocity(velocityLB);
-    para->setViscosity(viscosityLB);
+    para->setVelocityLB(velocityLB);
+    para->setViscosityLB(viscosityLB);
 
     para->setVelocityRatio(velocity/ velocityLB);
 
diff --git a/apps/gpu/LBM/gridGeneratorTest/gridGenerator.cpp b/apps/gpu/LBM/gridGeneratorTest/gridGenerator.cpp
index ecb43f56c79c25fc3aed6f404d717b8181fb0de1..745092f1711bd7cc523d1cd300fa785823634a3e 100644
--- a/apps/gpu/LBM/gridGeneratorTest/gridGenerator.cpp
+++ b/apps/gpu/LBM/gridGeneratorTest/gridGenerator.cpp
@@ -182,8 +182,8 @@ void multipleLevel(const std::string& configPath)
 
             para->setPrintFiles(true);
     
-            para->setVelocity( vx );
-            para->setViscosity( ( vx * D / dx ) / Re );
+            para->setVelocityLB( vx );
+            para->setViscosityLB( ( vx * D / dx ) / Re );
 
             para->setVelocityRatio(1.0);
 
@@ -270,8 +270,8 @@ void multipleLevel(const std::string& configPath)
 
             para->setPrintFiles(true);
     
-            para->setVelocity( vx );
-            para->setViscosity( ( vx * L / dx ) / Re );
+            para->setVelocityLB( vx );
+            para->setViscosityLB( ( vx * L / dx ) / Re );
 
             //para->setVelocityRatio(1.0 / velocityLB);
             para->setVelocityRatio(1.0);
@@ -376,8 +376,8 @@ void multipleLevel(const std::string& configPath)
 
             para->setPrintFiles(true);
     
-            para->setVelocity( vx );
-            para->setViscosity( ( vx * L / dx ) / Re );
+            para->setVelocityLB( vx );
+            para->setViscosityLB( ( vx * L / dx ) / Re );
 
             para->setVelocityRatio(1.0);
 
@@ -456,8 +456,8 @@ void multipleLevel(const std::string& configPath)
 
             para->setPrintFiles(true);
     
-            para->setVelocity( vx );
-            para->setViscosity( ( vx * L / dx ) / Re );
+            para->setVelocityLB( vx );
+            para->setViscosityLB( ( vx * L / dx ) / Re );
 
             para->setVelocityRatio(1.0);
 
@@ -653,8 +653,8 @@ void multipleLevel(const std::string& configPath)
 
             para->setPrintFiles(true);
     
-            para->setVelocity( vx );
-            para->setViscosity( ( vx * D / dx ) / Re );
+            para->setVelocityLB( vx );
+            para->setViscosityLB( ( vx * D / dx ) / Re );
 
             para->setVelocityRatio(1.0);
 
diff --git a/apps/gpu/LBM/lbmTest/main.cpp b/apps/gpu/LBM/lbmTest/main.cpp
index bd33f7c987934b73f41e2057da5e42655c406687..0fbd84b8291178988df8589dc32cad3a43aa10bf 100644
--- a/apps/gpu/LBM/lbmTest/main.cpp
+++ b/apps/gpu/LBM/lbmTest/main.cpp
@@ -106,8 +106,8 @@ void setParameters(std::shared_ptr<Parameter> para, std::unique_ptr<input::Input
     para->setTemperatureInit(StringUtil::toFloat(input->getValue("Temp")));
     para->setTemperatureBC(StringUtil::toFloat(input->getValue("TempBC")));
     //////////////////////////////////////////////////////////////////////////
-    para->setViscosity(StringUtil::toFloat(input->getValue("Viscosity_LB")));
-    para->setVelocity(StringUtil::toFloat(input->getValue("Velocity_LB")));
+    para->setViscosityLB(StringUtil::toFloat(input->getValue("Viscosity_LB")));
+    para->setVelocityLB(StringUtil::toFloat(input->getValue("Velocity_LB")));
     para->setViscosityRatio(StringUtil::toFloat(input->getValue("Viscosity_Ratio_World_to_LB")));
     para->setVelocityRatio(StringUtil::toFloat(input->getValue("Velocity_Ratio_World_to_LB")));
     para->setDensityRatio(StringUtil::toFloat(input->getValue("Density_Ratio_World_to_LB")));
diff --git a/apps/gpu/LBM/metisTest/main.cpp b/apps/gpu/LBM/metisTest/main.cpp
index 679119509d66737f3db132b676c73de5a146b779..cccf419062208071d7bd70cf449415acd009768d 100644
--- a/apps/gpu/LBM/metisTest/main.cpp
+++ b/apps/gpu/LBM/metisTest/main.cpp
@@ -106,8 +106,8 @@ void setParameters(std::shared_ptr<Parameter> para, std::unique_ptr<input::Input
     para->setTemperatureInit(StringUtil::toFloat(input->getValue("Temp")));
     para->setTemperatureBC(StringUtil::toFloat(input->getValue("TempBC")));
     //////////////////////////////////////////////////////////////////////////
-    para->setViscosity(StringUtil::toFloat(input->getValue("Viscosity_LB")));
-    para->setVelocity(StringUtil::toFloat(input->getValue("Velocity_LB")));
+    para->setViscosityLB(StringUtil::toFloat(input->getValue("Viscosity_LB")));
+    para->setVelocityLB(StringUtil::toFloat(input->getValue("Velocity_LB")));
     para->setViscosityRatio(StringUtil::toFloat(input->getValue("Viscosity_Ratio_World_to_LB")));
     para->setVelocityRatio(StringUtil::toFloat(input->getValue("Velocity_Ratio_World_to_LB")));
     para->setDensityRatio(StringUtil::toFloat(input->getValue("Density_Ratio_World_to_LB")));
diff --git a/apps/gpu/tests/NumericalTests/Utilities/VirtualFluidSimulationFactory/VirtualFluidSimulationFactoryImp.cpp b/apps/gpu/tests/NumericalTests/Utilities/VirtualFluidSimulationFactory/VirtualFluidSimulationFactoryImp.cpp
index 243521cf6a05899dfda957fd247e9cb6598d36a9..ca14490ae8e1656287d4a54e00bd72e781f5874b 100644
--- a/apps/gpu/tests/NumericalTests/Utilities/VirtualFluidSimulationFactory/VirtualFluidSimulationFactoryImp.cpp
+++ b/apps/gpu/tests/NumericalTests/Utilities/VirtualFluidSimulationFactory/VirtualFluidSimulationFactoryImp.cpp
@@ -42,8 +42,8 @@ std::shared_ptr<Parameter> VirtualFluidSimulationFactoryImp::makeParameter(std::
 	para->setTOut(simPara->getTimeStepLength());
 	para->setTStartOut(1);
 
-	para->setViscosity(simPara->getViscosity());
-	para->setVelocity(simPara->getMaxVelocity());
+	para->setViscosityLB(simPara->getViscosity());
+	para->setVelocityLB(simPara->getMaxVelocity());
 	para->setViscosityRatio(1.0);
 	para->setVelocityRatio(1.0);
 	para->setDensityRatio(1.0);
diff --git a/src/cuda/CudaGrid.cpp b/src/cuda/CudaGrid.cpp
index 9590452e107d17e69dd77b7159c20ca009f01a4c..5c514fcaddfd5b9a68c3ce20fd92e193700890d9 100644
--- a/src/cuda/CudaGrid.cpp
+++ b/src/cuda/CudaGrid.cpp
@@ -7,17 +7,7 @@ namespace vf::cuda
 
 CudaGrid::CudaGrid(unsigned int numberOfThreads, unsigned int numberOfEntities)
 {
-    unsigned int Grid = (numberOfEntities / numberOfThreads) + 1;
-    unsigned int Grid1, Grid2;
-    if (Grid > 512) {
-        Grid1 = 512;
-        Grid2 = (Grid / Grid1) + 1;
-    } else {
-        Grid1 = 1;
-        Grid2 = Grid;
-    }
-    
-    grid = dim3(Grid1, Grid2);
+    grid = getCudaGrid( numberOfThreads, numberOfEntities);
     threads = dim3(numberOfThreads, 1, 1);
 }
 
diff --git a/src/gpu/GridGenerator/geometries/BoundingBox/BoundingBoxTest.cpp b/src/gpu/GridGenerator/geometries/BoundingBox/BoundingBoxTest.cpp
index 67ac560b91f0ae712dfbc67d22676d8af8fc9afe..3d94d8fbcadb7d9803184885e67d1cb31e906d86 100644
--- a/src/gpu/GridGenerator/geometries/BoundingBox/BoundingBoxTest.cpp
+++ b/src/gpu/GridGenerator/geometries/BoundingBox/BoundingBoxTest.cpp
@@ -31,15 +31,15 @@ TEST(BoundingBoxExactTest, findMinMaxFromTriangle)
     Vertex normal = Vertex(0.0f, 0.0f, 0.0f);
     Triangle t = Triangle(v1, v2, v3, normal);
 
-	box.setMinMax(t);
-
-	EXPECT_THAT(box.minX, RealEq(minX));
-	EXPECT_THAT(box.minY, RealEq(minY));
-	EXPECT_THAT(box.minZ, RealEq(minZ));
-	
-	EXPECT_THAT(box.maxX, RealEq(maxX));
-	EXPECT_THAT(box.maxY, RealEq(maxY));
-	EXPECT_THAT(box.maxZ, RealEq(maxZ));
+    box.setMinMax(t);
+
+    EXPECT_THAT(box.minX, RealEq(minX));
+    EXPECT_THAT(box.minY, RealEq(minY));
+    EXPECT_THAT(box.minZ, RealEq(minZ));
+    
+    EXPECT_THAT(box.maxX, RealEq(maxX));
+    EXPECT_THAT(box.maxY, RealEq(maxY));
+    EXPECT_THAT(box.maxZ, RealEq(maxZ));
 }
 
 TEST(BoundingBoxTest, isInside_true)
diff --git a/src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.cpp b/src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.cpp
index a11384887074aa6b42bff77dd8b7ee1ade8fc9e0..883ca0deaf34f45e4608c4e59908b4562932db77 100644
--- a/src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.cpp
+++ b/src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.cpp
@@ -58,9 +58,9 @@ TriangularMesh* TriangularMesh::make(const std::string& fileName, const std::vec
 
 TriangularMesh::TriangularMesh(const std::string& input, const BoundingBox& box)
 {
-	this->triangleVec = STLReader::readSTL(box, input);
-	initalizeDataFromTriangles();
-	this->findNeighbors();
+    this->triangleVec = STLReader::readSTL(box, input);
+    initalizeDataFromTriangles();
+    this->findNeighbors();
 }
 
 TriangularMesh::TriangularMesh(const std::string& inputPath, const std::vector<uint> ignorePatches)
@@ -76,12 +76,7 @@ TriangularMesh::TriangularMesh(const std::string& inputPath, const std::vector<u
 
 TriangularMesh::TriangularMesh()
 {
-	this->minmax = BoundingBox::makeInvalidMinMaxBox();  // blame Lenz
-}
-
-TriangularMesh::~TriangularMesh()
-{
-
+    this->minmax = BoundingBox::makeInvalidMinMaxBox();  // blame Lenz
 }
 
 Object* TriangularMesh::clone() const
@@ -100,12 +95,12 @@ uint TriangularMesh::getNumberOfTriangles() const
 
 void TriangularMesh::findNeighbors()
 {
-	*logging::out << logging::Logger::INFO_INTERMEDIATE << "start finding neighbors ...\n";
+    *logging::out << logging::Logger::INFO_INTERMEDIATE << "start finding neighbors ...\n";
 
     auto t = Timer::makeStart();
 
-	TriangleNeighborFinder finder(triangles, size);
-	finder.fillWithNeighborAngles(this);
+    TriangleNeighborFinder finder(triangles, size);
+    finder.fillWithNeighborAngles(this);
 
     t->end();
 
@@ -114,19 +109,19 @@ void TriangularMesh::findNeighbors()
 
 void TriangularMesh::setTriangles(std::vector<Triangle> triangles)
 {
-	this->triangleVec = triangles;
-	initalizeDataFromTriangles();
+    this->triangleVec = triangles;
+    initalizeDataFromTriangles();
 }
 
 void TriangularMesh::setMinMax(BoundingBox minmax)
 {
-	this->minmax = minmax;
+    this->minmax = minmax;
 }
 
 void TriangularMesh::initalizeDataFromTriangles()
 {
-	this->triangles = triangleVec.data();
-	this->size = long(triangleVec.size());
+    this->triangles = triangleVec.data();
+    this->size = long(triangleVec.size());
 
     for (std::size_t i = 0; i < (size_t)this->size; i++) {
         this->minmax.setMinMax(this->triangleVec[i]);
@@ -201,7 +196,7 @@ void TriangularMesh::scale(double offset)
     auto averrageNormals = getAverrageNormalsPerVertex(trianglesPerVertex);
 
 
-    for (std::size_t vertexID = 0; vertexID < this->getNumberOfTriangles() * 3; vertexID++)
+    for (uint vertexID = 0; vertexID < this->getNumberOfTriangles() * 3; vertexID++)
     {
         int coordinatedID = finder.sortedToTriangles[vertexID][IDS::coordinateID];
         Vertex averrageNormal = averrageNormals[coordinatedID];
diff --git a/src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.h b/src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.h
index beb1c3a05b25904446e72a62196bcf6213fb0691..2e876e1d3c50b377ef6df9a8489fe8a189849594 100644
--- a/src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.h
+++ b/src/gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.h
@@ -33,7 +33,6 @@
 #ifndef TriangularMesh_h
 #define TriangularMesh_h
 
-#include <stdio.h>
 #include <vector>
 #include <string>
 #include <memory>
@@ -55,20 +54,20 @@ class TriangularMesh : public Object
 public:
 
     GRIDGENERATOR_EXPORT static TriangularMesh* make(const std::string& fileName, const std::vector<uint> ignorePatches = std::vector<uint>());
-	GRIDGENERATOR_EXPORT TriangularMesh();
+    GRIDGENERATOR_EXPORT TriangularMesh();
     GRIDGENERATOR_EXPORT TriangularMesh(const std::string& inputPath, const std::vector<uint> ignorePatches = std::vector<uint>());
-	GRIDGENERATOR_EXPORT TriangularMesh(const std::string& inputPath, const BoundingBox &box);
-	GRIDGENERATOR_EXPORT ~TriangularMesh();
+    GRIDGENERATOR_EXPORT TriangularMesh(const std::string& inputPath, const BoundingBox &box);
+    GRIDGENERATOR_EXPORT ~TriangularMesh() override = default;
 
     GRIDGENERATOR_EXPORT uint getNumberOfTriangles() const;
 
-	GRIDGENERATOR_EXPORT void setTriangles(std::vector<Triangle> triangles);
-	GRIDGENERATOR_EXPORT void setMinMax(BoundingBox minmax);
+    GRIDGENERATOR_EXPORT void setTriangles(std::vector<Triangle> triangles);
+    GRIDGENERATOR_EXPORT void setMinMax(BoundingBox minmax);
 
-	std::vector<Triangle> triangleVec;
-	Triangle *triangles;
-	long size;
-	BoundingBox minmax;
+    std::vector<Triangle> triangleVec;
+    Triangle *triangles = nullptr;
+    long size = 0;
+    BoundingBox minmax;
 
     SPtr<GbTriFaceMesh3D> VF_GbTriFaceMesh3D;
 
@@ -81,8 +80,8 @@ public:
     GRIDGENERATOR_EXPORT void generateGbTriFaceMesh3D();
 
 private:
-	
-	void initalizeDataFromTriangles();
+
+    void initalizeDataFromTriangles();
 
     static std::vector<Vertex> getAverrageNormalsPerVertex(std::vector<std::vector<Triangle> > trianglesPerVertex);
     static void eliminateTriangleswithIdenticialNormal(std::vector<Triangle> &triangles);
@@ -110,4 +109,3 @@ public:
 
 
 #endif
-
diff --git a/src/gpu/GridGenerator/grid/BoundaryConditions/Side.h b/src/gpu/GridGenerator/grid/BoundaryConditions/Side.h
index c9ffd40b0aa8fc2b8da8b4d85de60faea6927117..6df6bfccc9a39b80de3ac43d057a03945d035b34 100644
--- a/src/gpu/GridGenerator/grid/BoundaryConditions/Side.h
+++ b/src/gpu/GridGenerator/grid/BoundaryConditions/Side.h
@@ -36,7 +36,7 @@
 #include <string>
 #include <vector>
 
-#include "global.h"
+#include "gpu/GridGenerator/global.h"
 
 #define X_INDEX 0
 #define Y_INDEX 1
diff --git a/src/gpu/GridGenerator/grid/Field.cpp b/src/gpu/GridGenerator/grid/Field.cpp
index 86985af60e1ca25c247b586dbc2f356c665a8875..1777712107cc630d8b717233c11c224e278f787a 100644
--- a/src/gpu/GridGenerator/grid/Field.cpp
+++ b/src/gpu/GridGenerator/grid/Field.cpp
@@ -115,7 +115,7 @@ bool Field::isStopperCoarseUnderFine(uint index) const
 
 bool Field::isStopperSolid(uint index) const
 {
-	return field[index] == STOPPER_SOLID;
+    return field[index] == STOPPER_SOLID;
 }
 
 bool Field::isStopper(uint index) const
diff --git a/src/gpu/GridGenerator/grid/Field.h b/src/gpu/GridGenerator/grid/Field.h
index 08fff6da7c5a3f431138dc5039b4d234493ae4b8..db664d8820cce5eaeb01d08ef4d37f1ed7bd7064 100644
--- a/src/gpu/GridGenerator/grid/Field.h
+++ b/src/gpu/GridGenerator/grid/Field.h
@@ -51,20 +51,20 @@ public:
     bool is(uint index, char type) const;
     bool isCoarseToFineNode(uint index) const;
     bool isFineToCoarseNode(uint index) const;
-	bool isFluid(uint index) const;
-	bool isInvalidSolid(uint index) const;
+    bool isFluid(uint index) const;
+    bool isInvalidSolid(uint index) const;
     bool isQ(uint index) const;
     bool isBoundaryConditionNode(uint index) const;
     bool isInvalidCoarseUnderFine(uint index) const;
     bool isStopperOutOfGrid(uint index) const;
     bool isStopperCoarseUnderFine(uint index) const;
-	bool isStopperSolid(uint index) const;
-	bool isStopper(uint index) const;
+    bool isStopperSolid(uint index) const;
+    bool isStopper(uint index) const;
     bool isInvalidOutOfGrid(uint index) const;
 
     void setFieldEntry(uint index, char val);
-	void setFieldEntryToFluid(uint index);
-	void setFieldEntryToInvalidSolid(uint index);
+    void setFieldEntryToFluid(uint index);
+    void setFieldEntryToInvalidSolid(uint index);
     void setFieldEntryToStopperOutOfGrid(uint index);
     void setFieldEntryToStopperOutOfGridBoundary(uint index);
     void setFieldEntryToStopperCoarseUnderFine(uint index);
@@ -72,7 +72,7 @@ public:
     void setFieldEntryToInvalidOutOfGrid(uint index);
 
 private:
-    char *field;
+    char *field = nullptr;
     uint size;
 };
 
diff --git a/src/gpu/GridGenerator/grid/GridImp.cpp b/src/gpu/GridGenerator/grid/GridImp.cpp
index f6afafcd521245222c33972dcb46a9e9b2879826..030f94344817521a6059fcd498cb7fc457a7d862 100644
--- a/src/gpu/GridGenerator/grid/GridImp.cpp
+++ b/src/gpu/GridGenerator/grid/GridImp.cpp
@@ -32,8 +32,6 @@
 //=======================================================================================
 #include "GridImp.h"
 
-#include <stdio.h>
-#include <time.h>
 #include <iostream>
 #include <omp.h>
 #include <sstream>
@@ -1265,7 +1263,7 @@ void GridImp::mesh(Triangle &triangle)
                     continue;
 
                 const Vertex point(x, y, z);
-                const int value = triangle.isUnderFace(point);
+                const char value = triangle.isUnderFace(point);
                 //setDebugPoint(index, value);
 
                 if (value == Q_DEPRECATED)
diff --git a/src/gpu/GridGenerator/grid/GridInterface.h b/src/gpu/GridGenerator/grid/GridInterface.h
index b5f71317e7755a8b6bcfe3da084e0fc9155642f8..713d495d4386e0fe743357a803b84be02c061561 100644
--- a/src/gpu/GridGenerator/grid/GridInterface.h
+++ b/src/gpu/GridGenerator/grid/GridInterface.h
@@ -47,9 +47,9 @@ public:
     void GRIDGENERATOR_EXPORT findBoundaryGridInterfaceCF(const uint& indexOnCoarseGrid, GridImp* coarseGrid, GridImp* fineGrid);
 
 
-	void GRIDGENERATOR_EXPORT findInterfaceCF_GKS(const uint& indexOnCoarseGrid, GridImp* coarseGrid, GridImp* fineGrid);
+    void GRIDGENERATOR_EXPORT findInterfaceCF_GKS(const uint& indexOnCoarseGrid, GridImp* coarseGrid, GridImp* fineGrid);
 
-	void GRIDGENERATOR_EXPORT findInterfaceFC(const uint& indexOnCoarseGrid, GridImp* coarseGrid, GridImp* fineGrid);
+    void GRIDGENERATOR_EXPORT findInterfaceFC(const uint& indexOnCoarseGrid, GridImp* coarseGrid, GridImp* fineGrid);
     void GRIDGENERATOR_EXPORT findOverlapStopper(const uint& indexOnCoarseGrid, GridImp* coarseGrid, GridImp* fineGrid);
     
     void GRIDGENERATOR_EXPORT findInvalidBoundaryNodes(const uint& indexOnCoarseGrid, GridImp* coarseGrid);
@@ -66,7 +66,7 @@ public:
         uint *fine, *coarse;
         uint numberOfEntries = 0;
         uint *offset;
-    } fc, cf;
+    } fc{}, cf{};
 
 
 private:
diff --git a/src/gpu/GridGenerator/io/GridVTKWriter/GridVTKWriter.cpp b/src/gpu/GridGenerator/io/GridVTKWriter/GridVTKWriter.cpp
index 35b3197ff7c3f37eb33809cc9a909f0085d2dffc..7f818f3217e682f21c2b41c62070924243fcb3b0 100644
--- a/src/gpu/GridGenerator/io/GridVTKWriter/GridVTKWriter.cpp
+++ b/src/gpu/GridGenerator/io/GridVTKWriter/GridVTKWriter.cpp
@@ -59,7 +59,7 @@ void GridVTKWriter::writeSparseGridToVTK(SPtr<Grid> grid, const std::string& nam
     writeVtkFile(grid);
 }
 
-void GridVTKWriter::writeGridToVTKXML(SPtr<Grid> grid, const std::string& name, WRITING_FORMAT format)
+void GridVTKWriter::writeGridToVTKXML(SPtr<Grid> grid, const std::string& name)
 {
     
     const uint chunkSize = 20000000; 
@@ -87,11 +87,9 @@ void GridVTKWriter::writeGridToVTKXML(SPtr<Grid> grid, const std::string& name,
 
         *logging::out << logging::Logger::INFO_INTERMEDIATE << "Write Grid to XML VTK (*.vtu) output file : " + name + "_Part_" + std::to_string(part) + "\n";
 
-        nodedatanames.push_back("types");
-        nodedatanames.push_back("sparse_id");
-        nodedatanames.push_back("matrix_id");
-        nodedatanames.push_back("isSendNode");
-        nodedatanames.push_back("isReceiveNode");
+        nodedatanames.emplace_back("types");
+        nodedatanames.emplace_back("sparse_id");
+        nodedatanames.emplace_back("matrix_id");
 
         nodedata.resize(nodedatanames.size());
 
@@ -113,14 +111,12 @@ void GridVTKWriter::writeGridToVTKXML(SPtr<Grid> grid, const std::string& name,
                     grid->transIndexToCoords(index, x, y, z);
 
                     nodeNumbers(xIndex, yIndex, zIndex) = nr++;
-                    nodes.push_back(UbTupleFloat3(float(x), float(y), float(z)));
+                    nodes.emplace_back(UbTupleFloat3(float(x), float(y), float(z)));
 
                     const char type = grid->getFieldEntry(grid->transCoordToIndex(x, y, z));
                     nodedata[0].push_back(type);
                     nodedata[1].push_back(grid->getSparseIndex(index));
                     nodedata[2].push_back(index);
-                    nodedata[3].push_back(grid->isSendNode(index));
-                    nodedata[4].push_back(grid->isReceiveNode(index));
                 }
             }
         }
@@ -150,7 +146,7 @@ void GridVTKWriter::writeGridToVTKXML(SPtr<Grid> grid, const std::string& name,
                     {
                         Cell cell(x, y, z, grid->getDelta());
                         //if (grid->nodeInCellIs(cell, INVALID_OUT_OF_GRID) || grid->nodeInCellIs(cell, INVALID_COARSE_UNDER_FINE))
-                        //	continue;
+                        // continue;
 
                         cells.push_back(makeUbTuple(uint(SWB), uint(SEB), uint(NEB), uint(NWB), uint(SWT), uint(SET), uint(NET), uint(NWT)));
                     }
@@ -163,7 +159,7 @@ void GridVTKWriter::writeGridToVTKXML(SPtr<Grid> grid, const std::string& name,
 
 }
 
-void GridVTKWriter::writeInterpolationCellsToVTKXML(SPtr<Grid> grid, SPtr<Grid> gridCoarse, const std::string& name, WRITING_FORMAT format)
+void GridVTKWriter::writeInterpolationCellsToVTKXML(SPtr<Grid> grid, SPtr<Grid> gridCoarse, const std::string& name)
 {
     std::vector<char> nodeInterpolationCellType( grid->getSize() );
     for( auto& type : nodeInterpolationCellType ) type = -1;
@@ -202,79 +198,79 @@ void GridVTKWriter::writeInterpolationCellsToVTKXML(SPtr<Grid> grid, SPtr<Grid>
         }
     }
 
-	std::vector<UbTupleFloat3> nodes;
-	std::vector<UbTupleInt8> cells;
-	std::vector<std::string> celldatanames;
-	std::vector< std::vector<double> > celldata;
-
-	celldatanames.push_back("InterpolationCells");
-    celldatanames.push_back("Offset");
-
-	celldata.resize(celldatanames.size());
-
-	CbArray3D<int> nodeNumbers(grid->getNumberOfNodesX(), grid->getNumberOfNodesY(), grid->getNumberOfNodesZ(), -1);
-	int nr = 0;
-
-	for (uint xIndex = 0; xIndex < grid->getNumberOfNodesX(); xIndex++)
-	{
-		for (uint yIndex = 0; yIndex < grid->getNumberOfNodesY(); yIndex++)
-		{
-			for (uint zIndex = 0; zIndex < grid->getNumberOfNodesZ(); zIndex++)
-			{
-				real x, y, z;
-				uint index = 
-					  grid->getNumberOfNodesX() * grid->getNumberOfNodesY() * zIndex
-					+ grid->getNumberOfNodesX() *                             yIndex
-					+ xIndex;
-
-				grid->transIndexToCoords(index, x, y, z);
-
-				nodeNumbers(xIndex, yIndex, zIndex) = nr++;
-				nodes.push_back(UbTupleFloat3(float(x), float(y), float(z)));
-			}
-		}
-	}
-
-	int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
-	for (uint xIndex = 0; xIndex < grid->getNumberOfNodesX() - 1; xIndex++)
-	{
-		for (uint yIndex = 0; yIndex < grid->getNumberOfNodesY() - 1; yIndex++)
-		{
-			for (uint zIndex = 0; zIndex < grid->getNumberOfNodesZ() - 1; zIndex++)
-			{
-				real x, y, z;
-				uint index = grid->getNumberOfNodesX() * grid->getNumberOfNodesY() * zIndex
-					+ grid->getNumberOfNodesX() *                             yIndex
-					+ xIndex;
-
-				grid->transIndexToCoords(index, x, y, z);
-
-				if ((SWB = nodeNumbers(xIndex, yIndex, zIndex)) >= 0
-					&& (SEB = nodeNumbers(xIndex + 1, yIndex, zIndex)) >= 0
-					&& (NEB = nodeNumbers(xIndex + 1, yIndex + 1, zIndex)) >= 0
-					&& (NWB = nodeNumbers(xIndex, yIndex + 1, zIndex)) >= 0
-					&& (SWT = nodeNumbers(xIndex, yIndex, zIndex + 1)) >= 0
-					&& (SET = nodeNumbers(xIndex + 1, yIndex, zIndex + 1)) >= 0
-					&& (NET = nodeNumbers(xIndex + 1, yIndex + 1, zIndex + 1)) >= 0
-					&& (NWT = nodeNumbers(xIndex, yIndex + 1, zIndex + 1)) >= 0)
-				{
-					Cell cell(x, y, z, grid->getDelta());
-					//if (grid->nodeInCellIs(cell, INVALID_OUT_OF_GRID) || grid->nodeInCellIs(cell, INVALID_COARSE_UNDER_FINE))
-					//	continue;
-
-					cells.push_back(makeUbTuple(SWB, SEB, NEB, NWB, SWT, SET, NET, NWT));
-
-				    //const char type = grid->getFieldEntry(grid->transCoordToIndex(nodes[SWB].v1, nodes[SWB].v2.v1, nodes[SWB].v2.v2));
-				    //const char type = grid->getFieldEntry(grid->transCoordToIndex(val<1>(nodes[SWB]), val<2>(nodes[SWB]), val<3>(nodes[SWB])));
-				    const char type = nodeInterpolationCellType[ grid->transCoordToIndex(val<1>(nodes[SWB]), val<2>(nodes[SWB]), val<3>(nodes[SWB])) ];
+    std::vector<UbTupleFloat3> nodes;
+    std::vector<UbTupleInt8> cells;
+    std::vector<std::string> celldatanames;
+    std::vector< std::vector<double> > celldata;
+
+    celldatanames.emplace_back("InterpolationCells");
+    celldatanames.emplace_back("Offset");
+
+    celldata.resize(celldatanames.size());
+
+    CbArray3D<int> nodeNumbers(grid->getNumberOfNodesX(), grid->getNumberOfNodesY(), grid->getNumberOfNodesZ(), -1);
+    int nr = 0;
+
+    for (uint xIndex = 0; xIndex < grid->getNumberOfNodesX(); xIndex++)
+    {
+        for (uint yIndex = 0; yIndex < grid->getNumberOfNodesY(); yIndex++)
+        {
+            for (uint zIndex = 0; zIndex < grid->getNumberOfNodesZ(); zIndex++)
+            {
+                real x, y, z;
+                uint index = 
+                      grid->getNumberOfNodesX() * grid->getNumberOfNodesY() * zIndex
+                    + grid->getNumberOfNodesX() *                             yIndex
+                    + xIndex;
+
+                grid->transIndexToCoords(index, x, y, z);
+
+                nodeNumbers(xIndex, yIndex, zIndex) = nr++;
+                nodes.emplace_back(UbTupleFloat3(float(x), float(y), float(z)));
+            }
+        }
+    }
+
+    int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
+    for (uint xIndex = 0; xIndex < grid->getNumberOfNodesX() - 1; xIndex++)
+    {
+        for (uint yIndex = 0; yIndex < grid->getNumberOfNodesY() - 1; yIndex++)
+        {
+            for (uint zIndex = 0; zIndex < grid->getNumberOfNodesZ() - 1; zIndex++)
+            {
+                real x, y, z;
+                uint index = grid->getNumberOfNodesX() * grid->getNumberOfNodesY() * zIndex
+                    + grid->getNumberOfNodesX() *                             yIndex
+                    + xIndex;
+
+                grid->transIndexToCoords(index, x, y, z);
+
+                if ((SWB = nodeNumbers(xIndex, yIndex, zIndex)) >= 0
+                    && (SEB = nodeNumbers(xIndex + 1, yIndex, zIndex)) >= 0
+                    && (NEB = nodeNumbers(xIndex + 1, yIndex + 1, zIndex)) >= 0
+                    && (NWB = nodeNumbers(xIndex, yIndex + 1, zIndex)) >= 0
+                    && (SWT = nodeNumbers(xIndex, yIndex, zIndex + 1)) >= 0
+                    && (SET = nodeNumbers(xIndex + 1, yIndex, zIndex + 1)) >= 0
+                    && (NET = nodeNumbers(xIndex + 1, yIndex + 1, zIndex + 1)) >= 0
+                    && (NWT = nodeNumbers(xIndex, yIndex + 1, zIndex + 1)) >= 0)
+                {
+                    Cell cell(x, y, z, grid->getDelta());
+                    //if (grid->nodeInCellIs(cell, INVALID_OUT_OF_GRID) || grid->nodeInCellIs(cell, INVALID_COARSE_UNDER_FINE))
+                    //    continue;
+
+                    cells.push_back(makeUbTuple(SWB, SEB, NEB, NWB, SWT, SET, NET, NWT));
+
+                    //const char type = grid->getFieldEntry(grid->transCoordToIndex(nodes[SWB].v1, nodes[SWB].v2.v1, nodes[SWB].v2.v2));
+                    //const char type = grid->getFieldEntry(grid->transCoordToIndex(val<1>(nodes[SWB]), val<2>(nodes[SWB]), val<3>(nodes[SWB])));
+                    const char type = nodeInterpolationCellType[ grid->transCoordToIndex(val<1>(nodes[SWB]), val<2>(nodes[SWB]), val<3>(nodes[SWB])) ];
                     const char offset = nodeOffset               [ grid->transCoordToIndex(val<1>(nodes[SWB]), val<2>(nodes[SWB]), val<3>(nodes[SWB])) ];
 
                     celldata[0].push_back( type );
                     celldata[1].push_back( offset );
-				}
-			}
-		}
-	}
+                }
+            }
+        }
+    }
     WbWriterVtkXmlBinary::getInstance()->writeOctsWithCellData(name, nodes, cells, celldatanames, celldata);
 }
 
@@ -323,18 +319,18 @@ void GridVTKWriter::openFile(const std::string& name, const std::string& mode)
 void GridVTKWriter::closeFile()
 {
     GridVTKWriter::end_line();
-	fclose(file);
+    fclose(file);
 }
 
 void GridVTKWriter::writeHeader()
 {
-	fprintf(file, "# vtk DataFile Version 3.0\n");
-	fprintf(file, "by MeshGenerator\n");
-	if (isBinaryWritingFormat())
-		fprintf(file, "BINARY\n");
-	else
-		fprintf(file, "ASCII\n");
-	fprintf(file, "DATASET UNSTRUCTURED_GRID\n");
+    fprintf(file, "# vtk DataFile Version 3.0\n");
+    fprintf(file, "by MeshGenerator\n");
+    if (isBinaryWritingFormat())
+        fprintf(file, "BINARY\n");
+    else
+        fprintf(file, "ASCII\n");
+    fprintf(file, "DATASET UNSTRUCTURED_GRID\n");
 }
 
 void GridVTKWriter::writePoints(SPtr<Grid> grid)
@@ -360,34 +356,34 @@ void GridVTKWriter::writePoints(SPtr<Grid> grid)
 
 void GridVTKWriter::writeCells(const unsigned int &size)
 {
-	fprintf(file, "\nCELLS %d %d\n", size, size * 2);
-	for (unsigned int i = 0; i < size; ++i)
-	{
-		if (isBinaryWritingFormat()){
-			write_int(1);
-			write_int(i);
-		}
-		else
-			fprintf(file, "1 %d\n", i);
-	}
-
-	fprintf(file, "\nCELL_TYPES %d\n", size);
-	for (unsigned int i = 0; i < size; ++i)
-	{
-		if (isBinaryWritingFormat())
-			write_int(1);
-		else
-			fprintf(file, "1 ");
-	}
-	if (!isBinaryWritingFormat())
+    fprintf(file, "\nCELLS %u %u\n", size, size * 2);
+    for (unsigned int i = 0; i < size; ++i)
+    {
+        if (isBinaryWritingFormat()){
+            write_int(1);
+            write_int(i);
+        }
+        else
+            fprintf(file, "1 %u\n", i);
+    }
+
+    fprintf(file, "\nCELL_TYPES %u\n", size);
+    for (unsigned int i = 0; i < size; ++i)
+    {
+        if (isBinaryWritingFormat())
+            write_int(1);
+        else
+            fprintf(file, "1 ");
+    }
+    if (!isBinaryWritingFormat())
         GridVTKWriter::end_line();
 }
 
 void GridVTKWriter::writeTypeHeader(const unsigned int &size)
 {
-	fprintf(file, "\nPOINT_DATA %d\n", size);
-	fprintf(file, "SCALARS type int\n");
-	fprintf(file, "LOOKUP_TABLE default\n");
+    fprintf(file, "\nPOINT_DATA %u\n", size);
+    fprintf(file, "SCALARS type int\n");
+    fprintf(file, "LOOKUP_TABLE default\n");
 }
 
 void GridVTKWriter::writeTypes(SPtr<Grid> grid)
@@ -406,38 +402,38 @@ void GridVTKWriter::writeTypes(SPtr<Grid> grid)
 
 void GridVTKWriter::end_line()
 {
-	char str2[8] = "\n";
-	fprintf(file, "%s", str2);
+    char str2[8] = "\n";
+    fprintf(file, "%s", str2);
 }
 
 void GridVTKWriter::write_int(int val)
 {
-	force_big_endian((unsigned char *)&val);
-	fwrite(&val, sizeof(int), 1, file);
+    force_big_endian((unsigned char *)&val);
+    fwrite(&val, sizeof(int), 1, file);
 }
 
 void GridVTKWriter::write_float(float val)
 {
-	force_big_endian((unsigned char *)&val);
-	fwrite(&val, sizeof(float), 1, file);
+    force_big_endian((unsigned char *)&val);
+    fwrite(&val, sizeof(float), 1, file);
 }
 
 
 void GridVTKWriter::force_big_endian(unsigned char *bytes)
 {
-	bool shouldSwap = false;
-	int tmp1 = 1;
-	unsigned char *tmp2 = (unsigned char *)&tmp1;
-	if (*tmp2 != 0)
-		shouldSwap = true;
-
-	if (shouldSwap)
-	{
-		unsigned char tmp = bytes[0];
-		bytes[0] = bytes[3];
-		bytes[3] = tmp;
-		tmp = bytes[1];
-		bytes[1] = bytes[2];
-		bytes[2] = tmp;
-	}
+    bool shouldSwap = false;
+    int tmp1 = 1;
+    unsigned char *tmp2 = (unsigned char *)&tmp1;
+    if (*tmp2 != 0)
+        shouldSwap = true;
+
+    if (shouldSwap)
+    {
+        unsigned char tmp = bytes[0];
+        bytes[0] = bytes[3];
+        bytes[3] = tmp;
+        tmp = bytes[1];
+        bytes[1] = bytes[2];
+        bytes[2] = tmp;
+    }
 }
diff --git a/src/gpu/GridGenerator/io/GridVTKWriter/GridVTKWriter.h b/src/gpu/GridGenerator/io/GridVTKWriter/GridVTKWriter.h
index cf33df096a6e670b65f79831d59927e3d7cea389..f8458ae27bccd4b51db6b13544a708c4d8518051 100644
--- a/src/gpu/GridGenerator/io/GridVTKWriter/GridVTKWriter.h
+++ b/src/gpu/GridGenerator/io/GridVTKWriter/GridVTKWriter.h
@@ -46,12 +46,12 @@ class GRIDGENERATOR_EXPORT GridVTKWriter
 {
 public:
     static void writeSparseGridToVTK(SPtr<Grid> grid, const std::string& name, WRITING_FORMAT format = WRITING_FORMAT::ASCII);
-    static void writeGridToVTKXML(SPtr<Grid> grid, const std::string& name, WRITING_FORMAT format = WRITING_FORMAT::ASCII);
-    static void writeInterpolationCellsToVTKXML(SPtr<Grid> grid, SPtr<Grid> gridCoarse, const std::string& name, WRITING_FORMAT format = WRITING_FORMAT::ASCII);
+    static void writeGridToVTKXML(SPtr<Grid> grid, const std::string& name);
+    static void writeInterpolationCellsToVTKXML(SPtr<Grid> grid, SPtr<Grid> gridCoarse, const std::string& name);
 
 private:
-    GridVTKWriter() {}
-    ~GridVTKWriter() {}
+    GridVTKWriter() = default;
+    ~GridVTKWriter() = default;
 
     static FILE *file;
     static WRITING_FORMAT format;
diff --git a/src/gpu/GridGenerator/io/SimulationFileWriter/SimulationFileWriter.cpp b/src/gpu/GridGenerator/io/SimulationFileWriter/SimulationFileWriter.cpp
index 320a6e5fb7bb8e52a335722bca71d7d6a2a0c6de..23fb0f4e7f3e16702e9cb2459606986af1032e49 100644
--- a/src/gpu/GridGenerator/io/SimulationFileWriter/SimulationFileWriter.cpp
+++ b/src/gpu/GridGenerator/io/SimulationFileWriter/SimulationFileWriter.cpp
@@ -37,7 +37,6 @@
 #include <iomanip>
 #include <omp.h>
 #include <cmath>
-#include <stdint.h>
 
 #include "Core/Timer/Timer.h"
 
@@ -56,7 +55,7 @@ using namespace vf::gpu;
 /*#################################################################################*/
 /*---------------------------------public methods----------------------------------*/
 /*---------------------------------------------------------------------------------*/
-void SimulationFileWriter::write(std::string folder, SPtr<GridBuilder> builder, FILEFORMAT format)
+void SimulationFileWriter::write(const std::string& folder, SPtr<GridBuilder> builder, FILEFORMAT format)
 {
     SimulationFileWriter::folder = folder;
 
@@ -133,7 +132,7 @@ void SimulationFileWriter::openFiles(SPtr<GridBuilder> builder)
     qNames.push_back(path + simulationFileNames::bottomBoundaryQ);
     qNames.push_back(path + simulationFileNames::frontBoundaryQ);
     qNames.push_back(path + simulationFileNames::backBoundaryQ);
-	qNames.push_back(path + simulationFileNames::geomBoundaryQ);
+    qNames.push_back(path + simulationFileNames::geomBoundaryQ);
 
     std::vector<std::string> valueNames;
     valueNames.push_back(path + simulationFileNames::inletBoundaryValues);
@@ -142,7 +141,7 @@ void SimulationFileWriter::openFiles(SPtr<GridBuilder> builder)
     valueNames.push_back(path + simulationFileNames::bottomBoundaryValues);
     valueNames.push_back(path + simulationFileNames::frontBoundaryValues);
     valueNames.push_back(path + simulationFileNames::backBoundaryValues);
-	valueNames.push_back(path + simulationFileNames::geomBoundaryValues);
+    valueNames.push_back(path + simulationFileNames::geomBoundaryValues);
 
     for (int i = 0; i < QFILES; i++){
         SPtr<std::ofstream> outQ(new std::ofstream);
@@ -232,7 +231,7 @@ void SimulationFileWriter::writeLevelSize(uint numberOfNodes, FILEFORMAT format)
     const std::string zeroGeo = "16 ";
 
     if (format == FILEFORMAT::BINARY)
-	{
+    {
         //const uint zeroIndex = 0;
         //const uint zeroGeo   = 16;
 
@@ -258,7 +257,7 @@ void SimulationFileWriter::writeLevelSize(uint numberOfNodes, FILEFORMAT format)
         geoVecFile      << numberOfNodes << "\n" << zeroGeo  ;
     }
     else 
-	{
+    {
         xCoordFile      << numberOfNodes << "\n" << zeroIndex << "\n";
         yCoordFile      << numberOfNodes << "\n" << zeroIndex << "\n";
         zCoordFile      << numberOfNodes << "\n" << zeroIndex << "\n";
@@ -319,7 +318,7 @@ void SimulationFileWriter::writeCoordsNeighborsGeo(SPtr<GridBuilder> builder, in
     grid->transIndexToCoords(index, x, y, z);
 
     if (format == FILEFORMAT::BINARY)
-	{
+    {
         double tmpX = (double)x;
         double tmpY = (double)y;
         double tmpZ = (double)z;
@@ -342,7 +341,7 @@ void SimulationFileWriter::writeCoordsNeighborsGeo(SPtr<GridBuilder> builder, in
         geoVecFile.write((char*)&type, sizeof(unsigned int));
     }
     else 
-	{
+    {
         xCoordFile << x << "\n";
         yCoordFile << y << "\n";
         zCoordFile << z << "\n";
@@ -446,7 +445,7 @@ void SimulationFileWriter::writeGridInterfaceOffsetToFile(uint numberOfNodes, st
 std::vector<std::vector<std::vector<real> > > SimulationFileWriter::createBCVectors(SPtr<Grid> grid)
 {
     std::vector<std::vector<std::vector<real> > > qs;
-	qs.resize(QFILES);
+    qs.resize(QFILES);
     for (uint i = 0; i < grid->getSize(); i++)
     {
         real x, y, z;
@@ -472,9 +471,9 @@ void SimulationFileWriter::addShortQsToVector(int index, std::vector<std::vector
 
     for (int i = grid->getEndDirection(); i >= 0; i--)
     {
-		/*int qIndex = i * grid->getSize() + grid->getSparseIndex(index);
-		real q = grid->getDistribution()[qIndex];*/
-		real q = grid->getQValue(index, i);
+        /*int qIndex = i * grid->getSize() + grid->getSparseIndex(index);
+        real q = grid->getDistribution()[qIndex];*/
+        real q = grid->getQValue(index, i);
         if (q > 0) {
             //printf("Q%d (old:%d, new:%d), : %2.8f \n", i, coordsVec[index].matrixIndex, index, grid.d.f[i * grid.size + coordsVec[index].matrixIndex]);
             qKey += (uint32_t)pow(2, 26 - i);
@@ -499,10 +498,10 @@ void SimulationFileWriter::addQsToVector(int index, std::vector<std::vector<std:
     {
         //int qIndex = i * grid->getSize() + grid->getSparseIndex(index);
         //real q = grid->getDistribution()[qIndex];
-		real q = grid->getQValue(index, i);
+        real q = grid->getQValue(index, i);
         qNode.push_back(q);
-		if (q > 0)
-			printf("Q= %f; Index = %d \n", q, index);
+        if (q > 0)
+            printf("Q= %f; Index = %d \n", q, index);
             //qNode.push_back(q);
   //      else
   //          qNode.push_back(-1);
@@ -540,8 +539,8 @@ void SimulationFileWriter::writeBoundaryQsFile(SPtr<GridBuilder> builder)
   //  for (int rb = 0; rb < QFILES; rb++) {
   //      for (int index = 0; index < qFiles[rb].size(); index++) {
   //          //writeBoundary(qFiles[rb][index], rb);
-		//	writeBoundaryShort(qFiles[rb][index], rb);
-		//}
+        //    writeBoundaryShort(qFiles[rb][index], rb);
+        //}
   //  }
 
     SideType sides[] = {SideType::MX, SideType::PX, SideType::PZ, SideType::MZ, SideType::MY, SideType::PY, SideType::GEOMETRY};
@@ -593,18 +592,18 @@ void SimulationFileWriter::writeBoundary(std::vector<real> boundary, int rb)
 
 void SimulationFileWriter::writeBoundaryShort(std::vector<real> boundary, int rb)
 {
-	uint32_t key = *((uint32_t*)&boundary[boundary.size() - 2]);
-	int index = (int)boundary[boundary.size() - 1];
+    uint32_t key = *((uint32_t*)&boundary[boundary.size() - 2]);
+    int index = (int)boundary[boundary.size() - 1];
 
-	*qStreams[rb] << (index + 1) << " " << key;
+    *qStreams[rb] << (index + 1) << " " << key;
 
-	for (std::size_t i = 0; i < boundary.size() - 2; i++) {
-		*qStreams[rb] << " " << std::fixed << std::setprecision(16) << boundary[i];
-	}
-	*valueStreams[rb] << (index + 1) << " 0 0 0";
+    for (std::size_t i = 0; i < boundary.size() - 2; i++) {
+        *qStreams[rb] << " " << std::fixed << std::setprecision(16) << boundary[i];
+    }
+    *valueStreams[rb] << (index + 1) << " 0 0 0";
 
-	*qStreams[rb] << "\n";
-	*valueStreams[rb] << "\n";
+    *qStreams[rb] << "\n";
+    *valueStreams[rb] << "\n";
 }
 
 void SimulationFileWriter::writeBoundaryShort(SPtr<Grid> grid, SPtr<gg::BoundaryCondition> boundaryCondition, uint side)
diff --git a/src/gpu/GridGenerator/io/SimulationFileWriter/SimulationFileWriter.h b/src/gpu/GridGenerator/io/SimulationFileWriter/SimulationFileWriter.h
index 4a4552f74b69949865e233014d74ac7168b36b31..f3851abfd3372e5d3548cf7c0cd02344aa8acbaa 100644
--- a/src/gpu/GridGenerator/io/SimulationFileWriter/SimulationFileWriter.h
+++ b/src/gpu/GridGenerator/io/SimulationFileWriter/SimulationFileWriter.h
@@ -60,7 +60,7 @@ enum class FILEFORMAT
 class SimulationFileWriter : private NonCreatable
 {
 public:
-    GRIDGENERATOR_EXPORT static void write(std::string folder, SPtr<GridBuilder> builder, FILEFORMAT format);
+    GRIDGENERATOR_EXPORT static void write(const std::string& folder, SPtr<GridBuilder> builder, FILEFORMAT format);
 
 private:
     static void write(SPtr<GridBuilder> builder, FILEFORMAT format);
@@ -85,12 +85,12 @@ private:
     static void addQsToVector(int index, std::vector<std::vector<std::vector<real> > > &qs, SPtr<Grid> grid);
     static void fillRBForNode(int index, int direction, int directionSign, int rb, std::vector<std::vector<std::vector<real> > > &qs, SPtr<Grid> grid);
     static void writeBoundary(std::vector<real> boundary, int rb);
-	static void writeBoundaryShort(std::vector<real> boundary, int rb);
-	static void writeBoundaryShort(SPtr<Grid> grid, SPtr<gg::BoundaryCondition> boundaryCondition, uint side);
+    static void writeBoundaryShort(std::vector<real> boundary, int rb);
+    static void writeBoundaryShort(SPtr<Grid> grid, SPtr<gg::BoundaryCondition> boundaryCondition, uint side);
 
     static void writeCommunicationFiles(SPtr<GridBuilder> builder);
 
-	static void closeFiles();
+    static void closeFiles();
 
 
     static std::ofstream xCoordFile;
diff --git a/src/gpu/VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h b/src/gpu/VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h
index c2d374c7df97eb83d363417dc4a13b42b7312cab..095c2f4b16735eda966b58f0dc012386b617a3f5 100644
--- a/src/gpu/VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h
+++ b/src/gpu/VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactory.h
@@ -8,7 +8,7 @@
 
 #include "LBM/LB.h"
 #include "Parameter/Parameter.h"
-#include "grid/BoundaryConditions/Side.h"
+#include "gpu/GridGenerator/grid/BoundaryConditions/Side.h"
 
 struct LBMSimulationParameter;
 class Parameter;
@@ -29,7 +29,9 @@ public:
         VelocityCompressible,
         //! - VelocityAndPressureCompressible = interpolated velocity boundary condition, based on subgrid distances.
         //! Also sets the pressure to the bulk pressure. Can be combined with OutflowNonReflective
-        VelocityAndPressureCompressible
+        VelocityAndPressureCompressible,
+        //! - NotSpecified =  the user did not set a boundary condition
+        NotSpecified
     };
 
     //! \brief An enumeration for selecting a no-slip boundary condition
@@ -55,7 +57,9 @@ public:
         SlipCompressible,
         //! - SlipCompressible = interpolated slip boundary condition, based on subgrid distances.
         //! With turbulent viscosity -> para->setUseTurbulentViscosity(true) has to be set to true
-        SlipCompressibleTurbulentViscosity
+        SlipCompressibleTurbulentViscosity,
+        //! - NotSpecified =  the user did not set a boundary condition
+        NotSpecified
     };
 
     //! \brief An enumeration for selecting a pressure boundary condition
@@ -69,7 +73,9 @@ public:
         //! - PressureNonEquilibriumCompressible = pressure boundary condition based on non-equilibrium
         PressureNonEquilibriumCompressible,
         //! - OutflowNonReflective = outflow boundary condition, should be combined with VelocityAndPressureCompressible
-        OutflowNonReflective
+        OutflowNonReflective,
+        //! - NotSpecified =  the user did not set a boundary condition
+        NotSpecified
     };
 
     //! \brief An enumeration for selecting a stress boundary condition
@@ -77,7 +83,9 @@ public:
         //! - StressCompressible
         StressCompressible,
         //! - StressBounceBack
-        StressBounceBack
+        StressBounceBack,
+        //! - NotSpecified =  the user did not set a boundary condition
+        NotSpecified
     };
 
     // enum class OutflowBoundaryCondition {};  // TODO:
@@ -88,34 +96,35 @@ public:
     void setSlipBoundaryCondition(const BoundaryConditionFactory::SlipBC boundaryConditionType);
     void setPressureBoundaryCondition(const BoundaryConditionFactory::PressureBC boundaryConditionType);
     void setStressBoundaryCondition(const BoundaryConditionFactory::StressBC boundaryConditionType);
-    //!param boundaryConditionType: a velocity, no-slip or slip boundary condition
+    //! \brief set a boundary condition for the geometry
+    //! param boundaryConditionType: a velocity, no-slip or slip boundary condition
     //! \details suggestions for boundaryConditionType:
     //!
     //! - velocity: VelocityIncompressible, VelocityCompressible, VelocityAndPressureCompressible
     //!
-    //! - no-slip:  NoSlipBounceBack, NoSlipIncompressible, NoSlipCompressible, NoSlip3rdMomentsCompressible
+    //! - no-slip: NoSlipBounceBack, NoSlipIncompressible, NoSlipCompressible, NoSlip3rdMomentsCompressible
     //!
-    //! - slip:     SlipIncompressible
+    //! - slip: only use a slip boundary condition which sets the normals
     void setGeometryBoundaryCondition(const std::variant<VelocityBC, NoSlipBC, SlipBC> boundaryConditionType);
 
     // void setOutflowBoundaryCondition(...); // TODO:
     // https://git.rz.tu-bs.de/m.schoenherr/VirtualFluids_dev/-/issues/16
 
-    boundaryCondition getVelocityBoundaryConditionPost(bool isGeometryBC = false) const;
-    boundaryCondition getNoSlipBoundaryConditionPost(bool isGeometryBC = false) const;
-    boundaryCondition getSlipBoundaryConditionPost(bool isGeometryBC = false) const;
-    boundaryCondition getPressureBoundaryConditionPre() const;
-    boundaryCondition getGeometryBoundaryConditionPost() const;
+    [[nodiscard]] boundaryCondition getVelocityBoundaryConditionPost(bool isGeometryBC = false) const;
+    [[nodiscard]] boundaryCondition getNoSlipBoundaryConditionPost(bool isGeometryBC = false) const;
+    [[nodiscard]] boundaryCondition getSlipBoundaryConditionPost(bool isGeometryBC = false) const;
+    [[nodiscard]] boundaryCondition getPressureBoundaryConditionPre() const;
+    [[nodiscard]] boundaryCondition getGeometryBoundaryConditionPost() const;
 
     boundaryConditionPara getStressBoundaryConditionPost() const;
 
 private:
-    VelocityBC velocityBoundaryCondition;
+    VelocityBC velocityBoundaryCondition = VelocityBC::NotSpecified;
     NoSlipBC noSlipBoundaryCondition = NoSlipBC::NoSlipImplicitBounceBack;
-    SlipBC slipBoundaryCondition;
-    PressureBC pressureBoundaryCondition;
+    SlipBC slipBoundaryCondition = SlipBC::NotSpecified;
+    PressureBC pressureBoundaryCondition = PressureBC::NotSpecified;
     std::variant<VelocityBC, NoSlipBC, SlipBC> geometryBoundaryCondition  = NoSlipBC::NoSlipImplicitBounceBack;
-    StressBC stressBoundaryCondition;
+    StressBC stressBoundaryCondition = StressBC::NotSpecified;
 
 
     // OutflowBoundaryConditon outflowBC // TODO: https://git.rz.tu-bs.de/m.schoenherr/VirtualFluids_dev/-/issues/16
diff --git a/src/gpu/VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactoryTest.cpp b/src/gpu/VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactoryTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e3ba6624c3b43a49a0748666e27d1dde0a38ae7d
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/BoundaryConditions/BoundaryConditionFactoryTest.cpp
@@ -0,0 +1,51 @@
+#include <gmock/gmock-matchers.h>
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+#include "BoundaryConditionFactory.h"
+
+TEST(BoundaryConditionFactoryTest, defaultVelocityBC)
+{
+    auto bcFactory = BoundaryConditionFactory();
+    auto bc = bcFactory.getVelocityBoundaryConditionPost();
+    EXPECT_THAT(bc, testing::Eq(nullptr));
+    EXPECT_THROW(bc(nullptr, nullptr), std::bad_function_call);
+}
+
+TEST(BoundaryConditionFactoryTest, defaultNoSlipBC)
+{
+    auto bcFactory = BoundaryConditionFactory();
+    auto bc = bcFactory.getNoSlipBoundaryConditionPost();
+    EXPECT_NO_THROW(bc(nullptr, nullptr)); // empty lambda function should not throw
+}
+
+TEST(BoundaryConditionFactoryTest, defaultSlipBC)
+{
+    auto bcFactory = BoundaryConditionFactory();
+    auto bc = bcFactory.getSlipBoundaryConditionPost();
+    EXPECT_THAT(bc, testing::Eq(nullptr));
+    EXPECT_THROW(bc(nullptr, nullptr), std::bad_function_call);
+}
+
+TEST(BoundaryConditionFactoryTest, defaultPressureBC)
+{
+    auto bcFactory = BoundaryConditionFactory();
+    auto bc = bcFactory.getPressureBoundaryConditionPre();
+    EXPECT_THAT(bc, testing::Eq(nullptr));
+    EXPECT_THROW(bc(nullptr, nullptr), std::bad_function_call);
+}
+
+TEST(BoundaryConditionFactoryTest, defaultGeometryBC)
+{
+    auto bcFactory = BoundaryConditionFactory();
+    auto bc = bcFactory.getGeometryBoundaryConditionPost();
+    EXPECT_NO_THROW(bc(nullptr, nullptr)); // empty lambda function should not throw
+}
+
+TEST(BoundaryConditionFactoryTest, defaultStressBC)
+{
+    auto bcFactory = BoundaryConditionFactory();
+    auto bc = bcFactory.getStressBoundaryConditionPost();
+    EXPECT_THAT(bc, testing::Eq(nullptr));
+    EXPECT_THROW(bc(nullptr, nullptr, 0), std::bad_function_call);
+}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
index 5d036a4d94da7177e8ee76f93a27b355f562a50e..161813426f7a62ba4462237d1a6b5528516240f9 100644
--- a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
+++ b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
@@ -41,7 +41,7 @@ void UpdateGrid27::updateGrid(int level, unsigned int t)
 
     //////////////////////////////////////////////////////////////////////////
 
-    preCollisionBC(level, t);
+    this->preCollisionBC(level, t);
 
     //////////////////////////////////////////////////////////////////////////
     if( level != para->getFine() )
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.cpp b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.cpp
index da4ba9d8b4ad094437ec17b17e7ea653842bf83e..eea34e6fb7b53e16387d48077296d5df22a3e031 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.cpp
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.cpp
@@ -634,7 +634,7 @@ void GridReader::setNoSlipQs(std::shared_ptr<BoundaryQs> boundaryQ) const
 			this->printQSize("no slip", boundaryQ, level);
 			this->setSizeNoSlip(boundaryQ, level);
 			this->initalQStruct(para->getParH(level)->noSlipBC, boundaryQ, level);
-            cudaMemoryManager->cudaCopyWallBC(level);
+            cudaMemoryManager->cudaCopyNoSlipBC(level);
 		}
 	}
 }
@@ -766,7 +766,7 @@ void GridReader::setSizeNoSlip(std::shared_ptr<BoundaryQs> boundaryQ, unsigned i
 {
 	para->getParH(level)->noSlipBC.numberOfBCnodes = boundaryQ->getSize(level);
 	para->getParD(level)->noSlipBC.numberOfBCnodes = para->getParH(level)->noSlipBC.numberOfBCnodes;
-    cudaMemoryManager->cudaAllocWallBC(level);
+    cudaMemoryManager->cudaAllocNoSlipBC(level);
 }
 
 void GridReader::setSizeGeoQs(std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
index 705c50ecf38707fe817dfbd408035196da8a7c78..d29435659fd32e8f441028627c54ab50fb45ff25 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
@@ -17,16 +17,13 @@ using namespace vf::lbm::dir;
 
 GridGenerator::GridGenerator(std::shared_ptr<GridBuilder> builder, std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemoryManager> cudaMemoryManager, vf::gpu::Communicator& communicator)
 {
-	this->builder = builder;
+    this->builder = builder;
     this->para = para;
     this->cudaMemoryManager = cudaMemoryManager;
     this->indexRearrangement = std::make_unique<IndexRearrangementForStreams>(para, builder, communicator);
 }
 
-GridGenerator::~GridGenerator()
-{
-
-}
+GridGenerator::~GridGenerator() = default;
 
 void GridGenerator::initalGridInformations()
 {
@@ -46,22 +43,22 @@ void GridGenerator::initalGridInformations()
 void GridGenerator::allocArrays_CoordNeighborGeo()
 {
     const uint numberOfLevels = builder->getNumberOfGridLevels();
-	std::cout << "Number of Level: " << numberOfLevels << std::endl;
-	int numberOfNodesGlobal = 0;
-	std::cout << "Number of Nodes: " << std::endl;
-	
-	for (uint level = 0; level < numberOfLevels; level++) 
-	{
-		const int numberOfNodesPerLevel = builder->getNumberOfNodes(level) + 1;
-		numberOfNodesGlobal += numberOfNodesPerLevel;
-		std::cout << "Level " << level << " = " << numberOfNodesPerLevel << " Nodes" << std::endl;
-	
-		setNumberOfNodes(numberOfNodesPerLevel, level);
-	
-		cudaMemoryManager->cudaAllocCoord(level);
+    std::cout << "Number of Level: " << numberOfLevels << std::endl;
+    int numberOfNodesGlobal = 0;
+    std::cout << "Number of Nodes: " << std::endl;
+    
+    for (uint level = 0; level < numberOfLevels; level++) 
+    {
+        const int numberOfNodesPerLevel = builder->getNumberOfNodes(level) + 1;
+        numberOfNodesGlobal += numberOfNodesPerLevel;
+        std::cout << "Level " << level << " = " << numberOfNodesPerLevel << " Nodes" << std::endl;
+    
+        setNumberOfNodes(numberOfNodesPerLevel, level);
+    
+        cudaMemoryManager->cudaAllocCoord(level);
         cudaMemoryManager->cudaAllocSP(level);
         //cudaMemoryManager->cudaAllocF3SP(level);
-		cudaMemoryManager->cudaAllocNeighborWSB(level);
+        cudaMemoryManager->cudaAllocNeighborWSB(level);
 
         if(para->getUseTurbulentViscosity())
             cudaMemoryManager->cudaAllocTurbulentViscosity(level);
@@ -69,18 +66,18 @@ void GridGenerator::allocArrays_CoordNeighborGeo()
         if(para->getIsBodyForce())
             cudaMemoryManager->cudaAllocBodyForce(level);
 
-		builder->getNodeValues(
-			para->getParH(level)->coordinateX,
-			para->getParH(level)->coordinateY,
-			para->getParH(level)->coordinateZ,
-			para->getParH(level)->neighborX,
-			para->getParH(level)->neighborY,
-			para->getParH(level)->neighborZ,
-			para->getParH(level)->neighborInverse,
-			para->getParH(level)->typeOfGridNode,
-			level);
+        builder->getNodeValues(
+            para->getParH(level)->coordinateX,
+            para->getParH(level)->coordinateY,
+            para->getParH(level)->coordinateZ,
+            para->getParH(level)->neighborX,
+            para->getParH(level)->neighborY,
+            para->getParH(level)->neighborZ,
+            para->getParH(level)->neighborInverse,
+            para->getParH(level)->typeOfGridNode,
+            level);
 
-		setInitalNodeValues(numberOfNodesPerLevel, level);
+        setInitalNodeValues(numberOfNodesPerLevel, level);
 
         cudaMemoryManager->cudaCopyNeighborWSB(level);
         cudaMemoryManager->cudaCopySP(level);
@@ -89,9 +86,9 @@ void GridGenerator::allocArrays_CoordNeighborGeo()
             cudaMemoryManager->cudaCopyBodyForce(level);
 
         //std::cout << verifyNeighborIndices(level);
-	}
-	std::cout << "Number of Nodes: " << numberOfNodesGlobal << std::endl;
-	std::cout << "-----finish Coord, Neighbor, Geo------" << std::endl;
+    }
+    std::cout << "Number of Nodes: " << numberOfNodesGlobal << std::endl;
+    std::cout << "-----finish Coord, Neighbor, Geo------" << std::endl;
 }
 
 void GridGenerator::allocArrays_fluidNodeIndices() {
@@ -114,54 +111,57 @@ void GridGenerator::allocArrays_fluidNodeIndicesBorder() {
 
 void GridGenerator::allocArrays_BoundaryValues()
 {
-	std::cout << "------read BoundaryValues------" << std::endl;
+    std::cout << "------read BoundaryValues------" << std::endl;
+    int blocks = 0;
 
     for (uint level = 0; level < builder->getNumberOfGridLevels(); level++) {
         const auto numberOfPressureValues = int(builder->getPressureSize(level));
-
         std::cout << "size pressure level " << level << " : " << numberOfPressureValues << std::endl;
+
         ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-        para->getParH(level)->pressureBC.numberOfBCnodes = numberOfPressureValues;
-        para->getParD(level)->pressureBC.numberOfBCnodes = numberOfPressureValues;
-        para->getParH(level)->numberOfPressureBCnodesRead = numberOfPressureValues * para->getD3Qxx();
-        para->getParD(level)->numberOfPressureBCnodesRead = numberOfPressureValues * para->getD3Qxx();
+        para->getParH(level)->pressureBC.numberOfBCnodes = 0;
         if (numberOfPressureValues > 1)
         {
+            blocks = (numberOfPressureValues / para->getParH(level)->numberofthreads) + 1;
+            para->getParH(level)->pressureBC.numberOfBCnodes = blocks * para->getParH(level)->numberofthreads;
             cudaMemoryManager->cudaAllocPress(level);
             builder->getPressureValues(para->getParH(level)->pressureBC.RhoBC, para->getParH(level)->pressureBC.k, para->getParH(level)->pressureBC.kN, level);
             cudaMemoryManager->cudaCopyPress(level);
         }
+        para->getParD(level)->pressureBC.numberOfBCnodes = para->getParH(level)->pressureBC.numberOfBCnodes;
+        para->getParH(level)->numberOfPressureBCnodesRead = para->getParH(level)->pressureBC.numberOfBCnodes * para->getD3Qxx();
+        para->getParD(level)->numberOfPressureBCnodesRead = para->getParH(level)->pressureBC.numberOfBCnodes * para->getD3Qxx();
     }
 
     for (uint level = 0; level < builder->getNumberOfGridLevels(); level++) {
         const auto numberOfSlipValues = int(builder->getSlipSize(level));
-
         std::cout << "size slip level " << level << " : " << numberOfSlipValues << std::endl;
+
         ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-        para->getParH(level)->slipBC.numberOfBCnodes = numberOfSlipValues;
-        para->getParD(level)->slipBC.numberOfBCnodes = numberOfSlipValues;
-        para->getParH(level)->numberOfSlipBCnodesRead = numberOfSlipValues * para->getD3Qxx();
-        para->getParD(level)->numberOfSlipBCnodesRead = numberOfSlipValues * para->getD3Qxx();
+        para->getParH(level)->slipBC.numberOfBCnodes = 0;
         if (numberOfSlipValues > 1)
         {
+            blocks = (numberOfSlipValues / para->getParH(level)->numberofthreads) + 1;
+            para->getParH(level)->slipBC.numberOfBCnodes = blocks * para->getParH(level)->numberofthreads;
             cudaMemoryManager->cudaAllocSlipBC(level);
             builder->getSlipValues(para->getParH(level)->slipBC.normalX, para->getParH(level)->slipBC.normalY, para->getParH(level)->slipBC.normalZ, para->getParH(level)->slipBC.k, level);
             cudaMemoryManager->cudaCopySlipBC(level);
         }
+        para->getParD(level)->slipBC.numberOfBCnodes = para->getParH(level)->slipBC.numberOfBCnodes;
+        para->getParH(level)->numberOfSlipBCnodesRead = para->getParH(level)->slipBC.numberOfBCnodes * para->getD3Qxx();
+        para->getParD(level)->numberOfSlipBCnodesRead = para->getParH(level)->slipBC.numberOfBCnodes * para->getD3Qxx();
     }
 
     for (uint level = 0; level < builder->getNumberOfGridLevels(); level++) {
         const auto numberOfStressValues = int(builder->getStressSize(level));
-
         std::cout << "size stress level " << level << " : " << numberOfStressValues << std::endl;
-        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-        para->getParH(level)->stressBC.numberOfBCnodes = numberOfStressValues;
-        para->getParD(level)->stressBC.numberOfBCnodes = numberOfStressValues;
-        para->getParH(level)->numberOfStressBCnodesRead = numberOfStressValues * para->getD3Qxx();
-        para->getParD(level)->numberOfStressBCnodesRead = numberOfStressValues * para->getD3Qxx();
 
+        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        para->getParH(level)->stressBC.numberOfBCnodes = 0;
         if (numberOfStressValues > 1)
         {
+            blocks = (numberOfStressValues / para->getParH(level)->numberofthreads) + 1;
+            para->getParH(level)->stressBC.numberOfBCnodes = blocks * para->getParH(level)->numberofthreads;
             cudaMemoryManager->cudaAllocStressBC(level);
             cudaMemoryManager->cudaAllocWallModel(level, para->getHasWallModelMonitor());
             builder->getStressValues(   para->getParH(level)->stressBC.normalX,  para->getParH(level)->stressBC.normalY,  para->getParH(level)->stressBC.normalZ, 
@@ -174,6 +174,9 @@ void GridGenerator::allocArrays_BoundaryValues()
             cudaMemoryManager->cudaCopyStressBC(level);
             cudaMemoryManager->cudaCopyWallModel(level, para->getHasWallModelMonitor());
         }
+        para->getParD(level)->stressBC.numberOfBCnodes = para->getParH(level)->stressBC.numberOfBCnodes;
+        para->getParH(level)->numberOfStressBCnodesRead = para->getParH(level)->stressBC.numberOfBCnodes * para->getD3Qxx();
+        para->getParD(level)->numberOfStressBCnodesRead = para->getParH(level)->stressBC.numberOfBCnodes * para->getD3Qxx();
     }
     
 
@@ -181,17 +184,13 @@ void GridGenerator::allocArrays_BoundaryValues()
         const auto numberOfVelocityValues = int(builder->getVelocitySize(level));
         std::cout << "size velocity level " << level << " : " << numberOfVelocityValues << std::endl;
         ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-        int blocks = (numberOfVelocityValues / para->getParH(level)->numberofthreads) + 1;
-        para->getParH(level)->velocityBC.kArray = blocks * para->getParH(level)->numberofthreads;
-        para->getParD(level)->velocityBC.kArray = para->getParH(level)->velocityBC.kArray;
-        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-        para->getParH(level)->velocityBC.numberOfBCnodes = numberOfVelocityValues;
-        para->getParD(level)->velocityBC.numberOfBCnodes = numberOfVelocityValues;
-        para->getParH(level)->numberOfVeloBCnodesRead = numberOfVelocityValues * para->getD3Qxx();
-        para->getParD(level)->numberOfVeloBCnodesRead = numberOfVelocityValues * para->getD3Qxx();
+
+        para->getParH(level)->velocityBC.numberOfBCnodes = 0;
 
         if (numberOfVelocityValues > 1)
         {
+            blocks = (numberOfVelocityValues / para->getParH(level)->numberofthreads) + 1;
+            para->getParH(level)->velocityBC.numberOfBCnodes = blocks * para->getParH(level)->numberofthreads;
             ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
             cudaMemoryManager->cudaAllocVeloBC(level);
             ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -205,95 +204,85 @@ void GridGenerator::allocArrays_BoundaryValues()
             ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
             // advection - diffusion stuff
             if (para->getDiffOn()==true){
-            	//////////////////////////////////////////////////////////////////////////
-            	para->getParH(level)->TempVel.kTemp = numberOfVelocityValues;
-            	//cout << "Groesse kTemp = " << para->getParH(i)->TempPress.kTemp << endl;
-            	std::cout << "getTemperatureInit = " << para->getTemperatureInit() << std::endl;
-            	std::cout << "getTemperatureBC = " << para->getTemperatureBC() << std::endl;
-            	//////////////////////////////////////////////////////////////////////////
+                //////////////////////////////////////////////////////////////////////////
+                para->getParH(level)->TempVel.kTemp = para->getParH(level)->velocityBC.numberOfBCnodes;
+                //cout << "Groesse kTemp = " << para->getParH(i)->TempPress.kTemp << endl;
+                std::cout << "getTemperatureInit = " << para->getTemperatureInit() << std::endl;
+                std::cout << "getTemperatureBC = " << para->getTemperatureBC() << std::endl;
+                //////////////////////////////////////////////////////////////////////////
                 cudaMemoryManager->cudaAllocTempVeloBC(level);
-            	//cout << "nach alloc " << endl;
-            	//////////////////////////////////////////////////////////////////////////
-            	for (int m = 0; m < numberOfVelocityValues; m++)
-            	{
-            		para->getParH(level)->TempVel.temp[m]      = para->getTemperatureInit();
-            		para->getParH(level)->TempVel.tempPulse[m] = para->getTemperatureBC();
-            		para->getParH(level)->TempVel.velo[m]      = para->getVelocity();
-            		para->getParH(level)->TempVel.k[m]         = para->getParH(level)->velocityBC.k[m];
-            	}
-            	//////////////////////////////////////////////////////////////////////////
-            	//cout << "vor copy " << endl;
+                //cout << "nach alloc " << endl;
+                //////////////////////////////////////////////////////////////////////////
+                for (uint m = 0; m < para->getParH(level)->velocityBC.numberOfBCnodes; m++)
+                {
+                    para->getParH(level)->TempVel.temp[m]      = para->getTemperatureInit();
+                    para->getParH(level)->TempVel.tempPulse[m] = para->getTemperatureBC();
+                    para->getParH(level)->TempVel.velo[m]      = para->getVelocity();
+                    para->getParH(level)->TempVel.k[m]         = para->getParH(level)->velocityBC.k[m];
+                }
+                //////////////////////////////////////////////////////////////////////////
+                //cout << "vor copy " << endl;
                 cudaMemoryManager->cudaCopyTempVeloBCHD(level);
-            	//cout << "nach copy " << endl;
-            	//////////////////////////////////////////////////////////////////////////
+                //cout << "nach copy " << endl;
+                //////////////////////////////////////////////////////////////////////////
             }
             //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
         }
+        para->getParD(level)->velocityBC.numberOfBCnodes = para->getParH(level)->velocityBC.numberOfBCnodes;
+        para->getParH(level)->numberOfVeloBCnodesRead = para->getParH(level)->velocityBC.numberOfBCnodes * para->getD3Qxx();
+        para->getParD(level)->numberOfVeloBCnodesRead = para->getParH(level)->velocityBC.numberOfBCnodes * para->getD3Qxx();
     }
 
 
 
     if (builder->hasGeometryValues()) {
-        para->setGeometryValues(true);
-        for (uint i = 0; i < builder->getNumberOfGridLevels(); i++) {
-            int numberOfGeometryValues = builder->getGeometrySize(i);
-            std::cout << "size geometry values, Level " << i << " : " << numberOfGeometryValues << std::endl;
+        para->setUseGeometryValues(true);
+        for (uint level = 0; level < builder->getNumberOfGridLevels(); level++) {
+            int numberOfGeometryValues = builder->getGeometrySize(level);
+            std::cout << "size geometry values, Level " << level << " : " << numberOfGeometryValues << std::endl;
             ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-            para->getParH(i)->geometryBC.numberOfBCnodes = numberOfGeometryValues;
-            para->getParD(i)->geometryBC.numberOfBCnodes = numberOfGeometryValues;
+
+            para->getParH(level)->geometryBC.numberOfBCnodes = 0;
             if (numberOfGeometryValues > 0)
             {
-
+                blocks = (numberOfGeometryValues / para->getParH(level)->numberofthreads) + 1;
+                para->getParH(level)->geometryBC.numberOfBCnodes = blocks * para->getParH(level)->numberofthreads;
                 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-                cudaMemoryManager->cudaAllocGeomValuesBC(i);
+                cudaMemoryManager->cudaAllocGeomValuesBC(level);
                 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-                //Indexarray
 
-                builder->getGeometryValues(para->getParH(i)->geometryBC.Vx, para->getParH(i)->geometryBC.Vy, para->getParH(i)->geometryBC.Vz, i);
+                builder->getGeometryValues(para->getParH(level)->geometryBC.Vx, para->getParH(level)->geometryBC.Vy, para->getParH(level)->geometryBC.Vz, level);
 
                 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-                for (int m = 0; m < numberOfGeometryValues; m++)
+                for (uint m = 0; m < para->getParH(level)->geometryBC.numberOfBCnodes; m++)
                 {
-                    para->getParH(i)->geometryBC.Vx[m] = para->getParH(i)->geometryBC.Vx[m] / para->getVelocityRatio();
-                    para->getParH(i)->geometryBC.Vy[m] = para->getParH(i)->geometryBC.Vy[m] / para->getVelocityRatio();
-                    para->getParH(i)->geometryBC.Vz[m] = para->getParH(i)->geometryBC.Vz[m] / para->getVelocityRatio();
-                    //para->getParH(i)->geometryBC.Vx[m] = para->getParH(i)->geometryBC.Vx[m] / 100.0f;
-                    //para->getParH(i)->geometryBC.Vy[m] = para->getParH(i)->geometryBC.Vy[m] / 100.0f;
-                    //para->getParH(i)->geometryBC.Vz[m] = para->getParH(i)->geometryBC.Vz[m] / 100.0f;
-                    //para->getParH(i)->geometryBC.Vx[m] = 0.0f;
-                    //para->getParH(i)->geometryBC.Vy[m] = 0.0f;
-                    //para->getParH(i)->geometryBC.Vz[m] = 0.0f;
+                    para->getParH(level)->geometryBC.Vx[m] = para->getParH(level)->geometryBC.Vx[m] / para->getVelocityRatio();
+                    para->getParH(level)->geometryBC.Vy[m] = para->getParH(level)->geometryBC.Vy[m] / para->getVelocityRatio();
+                    para->getParH(level)->geometryBC.Vz[m] = para->getParH(level)->geometryBC.Vz[m] / para->getVelocityRatio();
                 }
-                //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-                ////T�st
-                //for (int m = 0; m < temp4; m++)
-                //{
-                //	para->getParH(i)->geometryBC.Vx[m] = para->getVelocity();//0.035f;
-                //	para->getParH(i)->geometryBC.Vy[m] = 0.0f;//para->getVelocity();//0.0f;
-                //	para->getParH(i)->geometryBC.Vz[m] = 0.0f;
-                //}
-                //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-                cudaMemoryManager->cudaCopyGeomValuesBC(i);
+                cudaMemoryManager->cudaCopyGeomValuesBC(level);
                 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
                 //// advection - diffusion stuff
                 //if (para->getDiffOn()==true){
-                //	//////////////////////////////////////////////////////////////////////////
-                //	para->getParH(i)->Temp.kTemp = temp4;
-                //	cout << "Groesse kTemp = " << para->getParH(i)->Temp.kTemp << std::endl;
-                //	//////////////////////////////////////////////////////////////////////////
-                //	para->cudaAllocTempNoSlipBC(i);
-                //	//////////////////////////////////////////////////////////////////////////
-                //	for (int m = 0; m < temp4; m++)
-                //	{
-                //		para->getParH(i)->Temp.temp[m] = para->getTemperatureInit();
-                //		para->getParH(i)->Temp.k[m]    = para->getParH(i)->geometryBC.k[m];
-                //	}
-                //	//////////////////////////////////////////////////////////////////////////
-                //	para->cudaCopyTempNoSlipBCHD(i);
-                //	//////////////////////////////////////////////////////////////////////////
+                //    //////////////////////////////////////////////////////////////////////////
+                //    para->getParH(i)->Temp.kTemp = temp4;
+                //    cout << "Groesse kTemp = " << para->getParH(i)->Temp.kTemp << std::endl;
+                //    //////////////////////////////////////////////////////////////////////////
+                //    para->cudaAllocTempNoSlipBC(i);
+                //    //////////////////////////////////////////////////////////////////////////
+                //    for (int m = 0; m < temp4; m++)
+                //    {
+                //        para->getParH(i)->Temp.temp[m] = para->getTemperatureInit();
+                //        para->getParH(i)->Temp.k[m]    = para->getParH(i)->geometryBC.k[m];
+                //    }
+                //    //////////////////////////////////////////////////////////////////////////
+                //    para->cudaCopyTempNoSlipBCHD(i);
+                //    //////////////////////////////////////////////////////////////////////////
                 //}
                 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
             }
+            para->getParD(level)->geometryBC.numberOfBCnodes = para->getParH(level)->geometryBC.numberOfBCnodes;
+
         }
     }//ende geo
 
@@ -731,11 +720,11 @@ void GridGenerator::initalValuesDomainDecompostion()
 
 void GridGenerator::allocArrays_BoundaryQs()
 {
-	std::cout << "------read BoundaryQs-------" << std::endl;
+    std::cout << "------read BoundaryQs-------" << std::endl;
 
 
     for (uint i = 0; i < builder->getNumberOfGridLevels(); i++) {
-        int numberOfPressureValues = (int)builder->getPressureSize(i);
+        const auto numberOfPressureValues = (int)builder->getPressureSize(i);
         if (numberOfPressureValues > 0)
         {
             std::cout << "size Pressure:  " << i << " : " << numberOfPressureValues << std::endl;
@@ -745,33 +734,7 @@ void GridGenerator::allocArrays_BoundaryQs()
             real* QQ = para->getParH(i)->pressureBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->pressureBC.numberOfBCnodes;
             QforBoundaryConditions Q;
-            Q.q27[E] = &QQ[E   *sizeQ];
-            Q.q27[W] = &QQ[W   *sizeQ];
-            Q.q27[N] = &QQ[N   *sizeQ];
-            Q.q27[S] = &QQ[S   *sizeQ];
-            Q.q27[T] = &QQ[T   *sizeQ];
-            Q.q27[B] = &QQ[B   *sizeQ];
-            Q.q27[NE] = &QQ[NE  *sizeQ];
-            Q.q27[SW] = &QQ[SW  *sizeQ];
-            Q.q27[SE] = &QQ[SE  *sizeQ];
-            Q.q27[NW] = &QQ[NW  *sizeQ];
-            Q.q27[TE] = &QQ[TE  *sizeQ];
-            Q.q27[BW] = &QQ[BW  *sizeQ];
-            Q.q27[BE] = &QQ[BE  *sizeQ];
-            Q.q27[TW] = &QQ[TW  *sizeQ];
-            Q.q27[TN] = &QQ[TN  *sizeQ];
-            Q.q27[BS] = &QQ[BS  *sizeQ];
-            Q.q27[BN] = &QQ[BN  *sizeQ];
-            Q.q27[TS] = &QQ[TS  *sizeQ];
-            Q.q27[REST] = &QQ[REST*sizeQ];
-            Q.q27[TNE] = &QQ[TNE *sizeQ];
-            Q.q27[TSW] = &QQ[TSW *sizeQ];
-            Q.q27[TSE] = &QQ[TSE *sizeQ];
-            Q.q27[TNW] = &QQ[TNW *sizeQ];
-            Q.q27[BNE] = &QQ[BNE *sizeQ];
-            Q.q27[BSW] = &QQ[BSW *sizeQ];
-            Q.q27[BSE] = &QQ[BSE *sizeQ];
-            Q.q27[BNW] = &QQ[BNW *sizeQ];
+            getPointersToBoundaryConditions(Q, QQ, sizeQ);
             
             builder->getPressureQs(Q.q27, i);
 
@@ -818,33 +781,7 @@ void GridGenerator::allocArrays_BoundaryQs()
             real* QQ = para->getParH(i)->slipBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->slipBC.numberOfBCnodes;
             QforBoundaryConditions Q;
-            Q.q27[E] = &QQ[E   *sizeQ];
-            Q.q27[W] = &QQ[W   *sizeQ];
-            Q.q27[N] = &QQ[N   *sizeQ];
-            Q.q27[S] = &QQ[S   *sizeQ];
-            Q.q27[T] = &QQ[T   *sizeQ];
-            Q.q27[B] = &QQ[B   *sizeQ];
-            Q.q27[NE] = &QQ[NE  *sizeQ];
-            Q.q27[SW] = &QQ[SW  *sizeQ];
-            Q.q27[SE] = &QQ[SE  *sizeQ];
-            Q.q27[NW] = &QQ[NW  *sizeQ];
-            Q.q27[TE] = &QQ[TE  *sizeQ];
-            Q.q27[BW] = &QQ[BW  *sizeQ];
-            Q.q27[BE] = &QQ[BE  *sizeQ];
-            Q.q27[TW] = &QQ[TW  *sizeQ];
-            Q.q27[TN] = &QQ[TN  *sizeQ];
-            Q.q27[BS] = &QQ[BS  *sizeQ];
-            Q.q27[BN] = &QQ[BN  *sizeQ];
-            Q.q27[TS] = &QQ[TS  *sizeQ];
-            Q.q27[REST] = &QQ[REST*sizeQ];
-            Q.q27[TNE] = &QQ[TNE *sizeQ];
-            Q.q27[TSW] = &QQ[TSW *sizeQ];
-            Q.q27[TSE] = &QQ[TSE *sizeQ];
-            Q.q27[TNW] = &QQ[TNW *sizeQ];
-            Q.q27[BNE] = &QQ[BNE *sizeQ];
-            Q.q27[BSW] = &QQ[BSW *sizeQ];
-            Q.q27[BSE] = &QQ[BSE *sizeQ];
-            Q.q27[BNW] = &QQ[BNW *sizeQ];
+            getPointersToBoundaryConditions(Q, QQ, sizeQ);
             
             builder->getSlipQs(Q.q27, i);
             ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -864,33 +801,7 @@ void GridGenerator::allocArrays_BoundaryQs()
             real* QQ = para->getParH(i)->stressBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->stressBC.numberOfBCnodes;
             QforBoundaryConditions Q;
-            Q.q27[E] = &QQ[E   *sizeQ];
-            Q.q27[W] = &QQ[W   *sizeQ];
-            Q.q27[N] = &QQ[N   *sizeQ];
-            Q.q27[S] = &QQ[S   *sizeQ];
-            Q.q27[T] = &QQ[T   *sizeQ];
-            Q.q27[B] = &QQ[B   *sizeQ];
-            Q.q27[NE] = &QQ[NE  *sizeQ];
-            Q.q27[SW] = &QQ[SW  *sizeQ];
-            Q.q27[SE] = &QQ[SE  *sizeQ];
-            Q.q27[NW] = &QQ[NW  *sizeQ];
-            Q.q27[TE] = &QQ[TE  *sizeQ];
-            Q.q27[BW] = &QQ[BW  *sizeQ];
-            Q.q27[BE] = &QQ[BE  *sizeQ];
-            Q.q27[TW] = &QQ[TW  *sizeQ];
-            Q.q27[TN] = &QQ[TN  *sizeQ];
-            Q.q27[BS] = &QQ[BS  *sizeQ];
-            Q.q27[BN] = &QQ[BN  *sizeQ];
-            Q.q27[TS] = &QQ[TS  *sizeQ];
-            Q.q27[REST] = &QQ[REST*sizeQ];
-            Q.q27[TNE] = &QQ[TNE *sizeQ];
-            Q.q27[TSW] = &QQ[TSW *sizeQ];
-            Q.q27[TSE] = &QQ[TSE *sizeQ];
-            Q.q27[TNW] = &QQ[TNW *sizeQ];
-            Q.q27[BNE] = &QQ[BNE *sizeQ];
-            Q.q27[BSW] = &QQ[BSW *sizeQ];
-            Q.q27[BSE] = &QQ[BSE *sizeQ];
-            Q.q27[BNW] = &QQ[BNW *sizeQ];
+            getPointersToBoundaryConditions(Q, QQ, sizeQ);
             
             builder->getStressQs(Q.q27, i);
             ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -910,34 +821,7 @@ void GridGenerator::allocArrays_BoundaryQs()
             real* QQ = para->getParH(i)->velocityBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->velocityBC.numberOfBCnodes;
             QforBoundaryConditions Q;
-            Q.q27[E] = &QQ[E   *sizeQ];
-            Q.q27[W] = &QQ[W   *sizeQ];
-            Q.q27[N] = &QQ[N   *sizeQ];
-            Q.q27[S] = &QQ[S   *sizeQ];
-            Q.q27[T] = &QQ[T   *sizeQ];
-            Q.q27[B] = &QQ[B   *sizeQ];
-            Q.q27[NE] = &QQ[NE  *sizeQ];
-            Q.q27[SW] = &QQ[SW  *sizeQ];
-            Q.q27[SE] = &QQ[SE  *sizeQ];
-            Q.q27[NW] = &QQ[NW  *sizeQ];
-            Q.q27[TE] = &QQ[TE  *sizeQ];
-            Q.q27[BW] = &QQ[BW  *sizeQ];
-            Q.q27[BE] = &QQ[BE  *sizeQ];
-            Q.q27[TW] = &QQ[TW  *sizeQ];
-            Q.q27[TN] = &QQ[TN  *sizeQ];
-            Q.q27[BS] = &QQ[BS  *sizeQ];
-            Q.q27[BN] = &QQ[BN  *sizeQ];
-            Q.q27[TS] = &QQ[TS  *sizeQ];
-            Q.q27[REST] = &QQ[REST*sizeQ];
-            Q.q27[TNE] = &QQ[TNE *sizeQ];
-            Q.q27[TSW] = &QQ[TSW *sizeQ];
-            Q.q27[TSE] = &QQ[TSE *sizeQ];
-            Q.q27[TNW] = &QQ[TNW *sizeQ];
-            Q.q27[BNE] = &QQ[BNE *sizeQ];
-            Q.q27[BSW] = &QQ[BSW *sizeQ];
-            Q.q27[BSE] = &QQ[BSE *sizeQ];
-            Q.q27[BNW] = &QQ[BNW *sizeQ];
-
+            getPointersToBoundaryConditions(Q, QQ, sizeQ);
             builder->getVelocityQs(Q.q27, i);
 
             if (para->getDiffOn()) {
@@ -993,37 +877,11 @@ void GridGenerator::allocArrays_BoundaryQs()
             real* QQ = para->getParH(i)->geometryBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->geometryBC.numberOfBCnodes;
             QforBoundaryConditions Q;
-            Q.q27[E] = &QQ[E   *sizeQ];
-            Q.q27[W] = &QQ[W   *sizeQ];
-            Q.q27[N] = &QQ[N   *sizeQ];
-            Q.q27[S] = &QQ[S   *sizeQ];
-            Q.q27[T] = &QQ[T   *sizeQ];
-            Q.q27[B] = &QQ[B   *sizeQ];
-            Q.q27[NE] = &QQ[NE  *sizeQ];
-            Q.q27[SW] = &QQ[SW  *sizeQ];
-            Q.q27[SE] = &QQ[SE  *sizeQ];
-            Q.q27[NW] = &QQ[NW  *sizeQ];
-            Q.q27[TE] = &QQ[TE  *sizeQ];
-            Q.q27[BW] = &QQ[BW  *sizeQ];
-            Q.q27[BE] = &QQ[BE  *sizeQ];
-            Q.q27[TW] = &QQ[TW  *sizeQ];
-            Q.q27[TN] = &QQ[TN  *sizeQ];
-            Q.q27[BS] = &QQ[BS  *sizeQ];
-            Q.q27[BN] = &QQ[BN  *sizeQ];
-            Q.q27[TS] = &QQ[TS  *sizeQ];
-            Q.q27[REST] = &QQ[REST*sizeQ];
-            Q.q27[TNE] = &QQ[TNE *sizeQ];
-            Q.q27[TSW] = &QQ[TSW *sizeQ];
-            Q.q27[TSE] = &QQ[TSE *sizeQ];
-            Q.q27[TNW] = &QQ[TNW *sizeQ];
-            Q.q27[BNE] = &QQ[BNE *sizeQ];
-            Q.q27[BSW] = &QQ[BSW *sizeQ];
-            Q.q27[BSE] = &QQ[BSE *sizeQ];
-            Q.q27[BNW] = &QQ[BNW *sizeQ];
+            getPointersToBoundaryConditions(Q, QQ, sizeQ);
             //////////////////////////////////////////////////////////////////
 
             builder->getGeometryQs(Q.q27, i);
-			//QDebugWriter::writeQValues(Q, para->getParH(i)->geometryBC.k, para->getParH(i)->geometryBC.numberOfBCnodes, "M:/TestGridGeneration/results/GeomGPU.dat");
+            //QDebugWriter::writeQValues(Q, para->getParH(i)->geometryBC.k, para->getParH(i)->geometryBC.numberOfBCnodes, "M:/TestGridGeneration/results/GeomGPU.dat");
             //////////////////////////////////////////////////////////////////
             for (int node_i = 0; node_i < numberOfGeometryNodes; node_i++)
             {
@@ -1031,10 +889,10 @@ void GridGenerator::allocArrays_BoundaryQs()
             }
             //for(int test = 0; test < 3; test++)
             //{
-            //	for (int tmp = 0; tmp < 27; tmp++)
-            //	{
-            //		cout <<"Kuhs: " << Q.q27[tmp][test]  << std::endl;
-            //	}
+            //    for (int tmp = 0; tmp < 27; tmp++)
+            //    {
+            //        cout <<"Kuhs: " << Q.q27[tmp][test]  << std::endl;
+            //    }
             //}
 
             ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -1063,7 +921,7 @@ void GridGenerator::allocArrays_BoundaryQs()
     }
 
 
-	std::cout << "-----finish BoundaryQs------" << std::endl;
+    std::cout << "-----finish BoundaryQs------" << std::endl;
 }
 
 void GridGenerator::allocArrays_OffsetScale()
@@ -1098,10 +956,10 @@ void GridGenerator::allocArrays_OffsetScale()
         para->getParD(level)->mem_size_kFC_off = sizeof(real)* para->getParD(level)->K_FC;
         ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
         //alloc
-		cudaMemoryManager->cudaAllocInterfaceCF(level);
-		cudaMemoryManager->cudaAllocInterfaceFC(level);
-		cudaMemoryManager->cudaAllocInterfaceOffCF(level);
-		cudaMemoryManager->cudaAllocInterfaceOffFC(level);
+        cudaMemoryManager->cudaAllocInterfaceCF(level);
+        cudaMemoryManager->cudaAllocInterfaceFC(level);
+        cudaMemoryManager->cudaAllocInterfaceOffCF(level);
+        cudaMemoryManager->cudaAllocInterfaceOffFC(level);
         ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
         //init
         builder->getOffsetCF(para->getParH(level)->offCF.xOffCF, para->getParH(level)->offCF.yOffCF, para->getParH(level)->offCF.zOffCF, level);
@@ -1116,49 +974,49 @@ void GridGenerator::allocArrays_OffsetScale()
         }
         ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
         //copy
-		cudaMemoryManager->cudaCopyInterfaceCF(level);
-		cudaMemoryManager->cudaCopyInterfaceFC(level);
-		cudaMemoryManager->cudaCopyInterfaceOffCF(level);
-		cudaMemoryManager->cudaCopyInterfaceOffFC(level);
+        cudaMemoryManager->cudaCopyInterfaceCF(level);
+        cudaMemoryManager->cudaCopyInterfaceFC(level);
+        cudaMemoryManager->cudaCopyInterfaceOffCF(level);
+        cudaMemoryManager->cudaCopyInterfaceOffFC(level);
         ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     }
 }
 
 void GridGenerator::setDimensions()
 {
-	//std::vector<int> localGridNX(1);
-	//std::vector<int> localGridNY(1);
-	//std::vector<int> localGridNZ(1);
+    //std::vector<int> localGridNX(1);
+    //std::vector<int> localGridNY(1);
+    //std::vector<int> localGridNZ(1);
 
-	//builder->getDimensions(localGridNX[0], localGridNY[0], localGridNZ[0], 0);
+    //builder->getDimensions(localGridNX[0], localGridNY[0], localGridNZ[0], 0);
 
-	//para->setGridX(localGridNX);
-	//para->setGridY(localGridNY);
-	//para->setGridZ(localGridNZ);
+    //para->setGridX(localGridNX);
+    //para->setGridY(localGridNY);
+    //para->setGridZ(localGridNZ);
 }
 
 void GridGenerator::setBoundingBox()
 {
-	std::vector<int> localGridNX(1);
-	std::vector<int> localGridNY(1);
-	std::vector<int> localGridNZ(1);
-	builder->getDimensions(localGridNX[0], localGridNY[0], localGridNZ[0], 0);
-
-	std::vector<real> minX, maxX, minY, maxY, minZ, maxZ;
-	minX.push_back(0);
-	minY.push_back(0);
-	minZ.push_back(0);
-
-	maxX.push_back((real)localGridNX[0]);
-	maxY.push_back((real)localGridNY[0]);
-	maxZ.push_back((real)localGridNZ[0]);
-
-	para->setMinCoordX(minX);
-	para->setMinCoordY(minY);
-	para->setMinCoordZ(minZ);
-	para->setMaxCoordX(maxX);
-	para->setMaxCoordY(maxY);
-	para->setMaxCoordZ(maxZ);
+    std::vector<int> localGridNX(1);
+    std::vector<int> localGridNY(1);
+    std::vector<int> localGridNZ(1);
+    builder->getDimensions(localGridNX[0], localGridNY[0], localGridNZ[0], 0);
+
+    std::vector<real> minX, maxX, minY, maxY, minZ, maxZ;
+    minX.push_back(0);
+    minY.push_back(0);
+    minZ.push_back(0);
+
+    maxX.push_back((real)localGridNX[0]);
+    maxY.push_back((real)localGridNY[0]);
+    maxZ.push_back((real)localGridNZ[0]);
+
+    para->setMinCoordX(minX);
+    para->setMinCoordY(minY);
+    para->setMinCoordZ(minZ);
+    para->setMaxCoordX(maxX);
+    para->setMaxCoordY(maxY);
+    para->setMaxCoordZ(maxZ);
 }
 
 void GridGenerator::initPeriodicNeigh(std::vector<std::vector<std::vector<uint> > > periodV, std::vector<std::vector<uint> > periodIndex, std::string way)
@@ -1258,3 +1116,33 @@ std::string GridGenerator::checkNeighbor(int level, real x, real y, real z, int
     }
     return oss.str();
 }
+
+void GridGenerator::getPointersToBoundaryConditions(QforBoundaryConditions& boundaryConditionStruct, real* subgridDistances, const unsigned int numberOfBCnodes){
+    boundaryConditionStruct.q27[E] =    &subgridDistances[E   * numberOfBCnodes];
+    boundaryConditionStruct.q27[W] =    &subgridDistances[W   * numberOfBCnodes];
+    boundaryConditionStruct.q27[N] =    &subgridDistances[N   * numberOfBCnodes];
+    boundaryConditionStruct.q27[S] =    &subgridDistances[S   * numberOfBCnodes];
+    boundaryConditionStruct.q27[T] =    &subgridDistances[T   * numberOfBCnodes];
+    boundaryConditionStruct.q27[B] =    &subgridDistances[B   * numberOfBCnodes];
+    boundaryConditionStruct.q27[NE] =   &subgridDistances[NE  * numberOfBCnodes];
+    boundaryConditionStruct.q27[SW] =   &subgridDistances[SW  * numberOfBCnodes];
+    boundaryConditionStruct.q27[SE] =   &subgridDistances[SE  * numberOfBCnodes];
+    boundaryConditionStruct.q27[NW] =   &subgridDistances[NW  * numberOfBCnodes];
+    boundaryConditionStruct.q27[TE] =   &subgridDistances[TE  * numberOfBCnodes];
+    boundaryConditionStruct.q27[BW] =   &subgridDistances[BW  * numberOfBCnodes];
+    boundaryConditionStruct.q27[BE] =   &subgridDistances[BE  * numberOfBCnodes];
+    boundaryConditionStruct.q27[TW] =   &subgridDistances[TW  * numberOfBCnodes];
+    boundaryConditionStruct.q27[TN] =   &subgridDistances[TN  * numberOfBCnodes];
+    boundaryConditionStruct.q27[BS] =   &subgridDistances[BS  * numberOfBCnodes];
+    boundaryConditionStruct.q27[BN] =   &subgridDistances[BN  * numberOfBCnodes];
+    boundaryConditionStruct.q27[TS] =   &subgridDistances[TS  * numberOfBCnodes];
+    boundaryConditionStruct.q27[REST] = &subgridDistances[REST* numberOfBCnodes];
+    boundaryConditionStruct.q27[TNE] =  &subgridDistances[TNE * numberOfBCnodes];
+    boundaryConditionStruct.q27[TSW] =  &subgridDistances[TSW * numberOfBCnodes];
+    boundaryConditionStruct.q27[TSE] =  &subgridDistances[TSE * numberOfBCnodes];
+    boundaryConditionStruct.q27[TNW] =  &subgridDistances[TNW * numberOfBCnodes];
+    boundaryConditionStruct.q27[BNE] =  &subgridDistances[BNE * numberOfBCnodes];
+    boundaryConditionStruct.q27[BSW] =  &subgridDistances[BSW * numberOfBCnodes];
+    boundaryConditionStruct.q27[BSE] =  &subgridDistances[BSE * numberOfBCnodes];
+    boundaryConditionStruct.q27[BNW] =  &subgridDistances[BNW * numberOfBCnodes];
+}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.h b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.h
index 991c59ca63dd77da974c044db363de16ebf6ff23..73cb0597067c63ba1ec40285dd9bede2349d5aa1 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.h
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.h
@@ -1,3 +1,35 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __         
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |        
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |        
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |        
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____    
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|   
+//      \    \  |    |   ________________________________________________________________    
+//       \    \ |    |  |  ______________________________________________________________|   
+//        \    \|    |  |  |         __          __     __     __     ______      _______    
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)   
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______    
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/   
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can 
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of 
+//  the License, or (at your option) any later version.
+//  
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT 
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
+//  for more details.
+//  
+//  You should have received a copy of the GNU General Public License along
+//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file GridGenerator.h
+//! \ingroup DataStructureInitializer
+//! \author Martin Schoenherr
+//=======================================================================================
 #ifndef GridReaderGenerator_H
 #define GridReaderGenerator_H
 
@@ -13,56 +45,62 @@ class Parameter;
 class GridBuilder;
 class IndexRearrangementForStreams;
 
+//! \class GridGenerator derived class of GridProvider
+//! \brief mapping the grid of grid generator to data structure for simulation
 class GridGenerator
-	: public GridProvider
+    : public GridProvider
 {
 private:
-	std::vector<std::string> channelDirections;
-	std::vector<std::string> channelBoundaryConditions;
+    //! \brief string vector with channel direction
+    std::vector<std::string> channelDirections;
+    //! \brief string vector with channel direction (boundary conditions)
+    std::vector<std::string> channelBoundaryConditions;
 
-	std::shared_ptr<GridBuilder> builder;
+    std::shared_ptr<GridBuilder> builder;
     std::unique_ptr<IndexRearrangementForStreams> indexRearrangement;
 
 public:
     VIRTUALFLUIDS_GPU_EXPORT GridGenerator(std::shared_ptr<GridBuilder> builder, std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemoryManager> cudaMemoryManager, vf::gpu::Communicator& communicator);
-	VIRTUALFLUIDS_GPU_EXPORT virtual ~GridGenerator();
+    VIRTUALFLUIDS_GPU_EXPORT ~GridGenerator() override;
 
-	void allocArrays_CoordNeighborGeo() override;
+    //! \brief allocates and initialized the data structures for Coordinates and node types
+    void allocArrays_CoordNeighborGeo() override;
+    //! \brief allocates and initialized the values at the boundary conditions
     void allocArrays_BoundaryValues() override;
-
-	void allocArrays_BoundaryQs() override;
+    //! \brief allocates and initialized the sub-grid distances at the boundary conditions
+    void allocArrays_BoundaryQs() override;
     void allocArrays_OffsetScale() override;
     void allocArrays_fluidNodeIndices() override;
     void allocArrays_fluidNodeIndicesBorder() override;
 
-	virtual void setDimensions() override;
-	virtual void setBoundingBox() override;
+    virtual void setDimensions() override;
+    virtual void setBoundingBox() override;
 
-	virtual void initPeriodicNeigh(std::vector<std::vector<std::vector<unsigned int> > > periodV, std::vector<std::vector<unsigned int> > periodIndex, std::string way) override;
-	
+    virtual void initPeriodicNeigh(std::vector<std::vector<std::vector<unsigned int> > > periodV, std::vector<std::vector<unsigned int> > periodIndex, std::string way) override;
+    
 private:
-	void setPressureValues(int channelSide) const;
-	void setPressRhoBC(int sizePerLevel, int level, int channelSide) const;
-
-	void setVelocityValues(int channelSide) const;
-	void setVelocity(int level, int sizePerLevel, int channelSide) const;
-
-	void setOutflowValues(int channelSide) const;
-	void setOutflow(int level, int sizePerLevel, int channelSide) const;
-
-	void setPressQs(int channelSide) const;
-	void setVelocityQs(int channelSide) const;
-	void setOutflowQs(int channelSide) const;
-	void setNoSlipQs(int channelSide) const;
-	void setGeoQs() const;
-	void modifyQElement(int channelSide, unsigned int level) const;
-
-	void initalQStruct(QforBoundaryConditions& Q,int channelSide, unsigned int level) const;
-	void printQSize(std::string bc,int channelSide, unsigned int level) const;
-	void setSizeNoSlip(int channelSide, unsigned int level) const;
-	void setSizeGeoQs(unsigned int level) const;
-	void setQ27Size(QforBoundaryConditions &Q, real* QQ, unsigned int sizeQ) const;
-	bool hasQs(int channelSide, unsigned int level) const;
+    void setPressureValues(int channelSide) const;
+    void setPressRhoBC(int sizePerLevel, int level, int channelSide) const;
+
+    void setVelocityValues(int channelSide) const;
+    void setVelocity(int level, int sizePerLevel, int channelSide) const;
+
+    void setOutflowValues(int channelSide) const;
+    void setOutflow(int level, int sizePerLevel, int channelSide) const;
+
+    void setPressQs(int channelSide) const;
+    void setVelocityQs(int channelSide) const;
+    void setOutflowQs(int channelSide) const;
+    void setNoSlipQs(int channelSide) const;
+    void setGeoQs() const;
+    void modifyQElement(int channelSide, unsigned int level) const;
+
+    void initalQStruct(QforBoundaryConditions& Q,int channelSide, unsigned int level) const;
+    void printQSize(std::string bc,int channelSide, unsigned int level) const;
+    void setSizeNoSlip(int channelSide, unsigned int level) const;
+    void setSizeGeoQs(unsigned int level) const;
+    void setQ27Size(QforBoundaryConditions &Q, real* QQ, unsigned int sizeQ) const;
+    bool hasQs(int channelSide, unsigned int level) const;
     
     void initalValuesDomainDecompostion();
 public:
@@ -70,10 +108,26 @@ public:
 
 
 private:
+    //! \brief verifies if there are invalid nodes, stopper nodes or wrong neighbors
     std::string verifyNeighborIndices(int level) const;
-    std::string verifyNeighborIndex(int level, int index, int &invalidNodes, int &stopperNodes, int &wrongNeighbors) const;
+    //! \brief verifies single neighbor index
+    //! \param index type integer
+    //! \param invalidNodes reference to invalid nodes
+    //! \param stopperNodes reference to stopper nodes
+    //! \param wrongNeighbors reference to wrong neighbors
+    std::string verifyNeighborIndex(int level, int index , int &invalidNodes, int &stopperNodes, int &wrongNeighbors) const;
+    //! \brief check the neighbors
+    //! \param x,y,z lattice node position
+    //! \param numberOfWrongNeighbors reference to the number of wrong neighbors
+    //! \param neighborIndex index of neighbor node
+    //! \param neighborX,neighborY,neighborZ neighbor lattice node position
+    //! \param direction type string
     std::string checkNeighbor(int level, real x, real y, real z, int index, int& numberOfWrongNeihgbors, int neighborIndex, real neighborX, real neighborY, real neighborZ, std::string direction) const;
-
+    //! \brief create the pointers in the struct for the BoundaryConditions from the boundary condition array
+    //! \param boundaryConditionStruct is a struct containing information about the boundary condition
+    //! \param subgridDistances is a pointer to an array containing the subgrid distances
+    //! \param numberOfBCnodes is the number of lattice nodes in the boundary condition
+    static void getPointersToBoundaryConditions(QforBoundaryConditions& boundaryConditionStruct, real* subgridDistances, const unsigned int numberOfBCnodes);
 };
 
 #endif
diff --git a/src/gpu/VirtualFluids_GPU/FindQ/DefineBCs.cpp b/src/gpu/VirtualFluids_GPU/FindQ/DefineBCs.cpp
index a67f1d987cb9636ee447c5f5acd1410c44cb6a62..a4e3123873ea82e0904efa31c9fc018dcc79c227 100644
--- a/src/gpu/VirtualFluids_GPU/FindQ/DefineBCs.cpp
+++ b/src/gpu/VirtualFluids_GPU/FindQ/DefineBCs.cpp
@@ -33,13 +33,13 @@ void findQ27(Parameter* para, CudaMemoryManager* cudaMemoryManager)
 	   para->getParD(lev)->noSlipBC.numberOfBCnodes = para->getParH(lev)->noSlipBC.numberOfBCnodes;
       printf("numberOfBCnodes= %d\n", para->getParH(lev)->noSlipBC.numberOfBCnodes);
 
-	  cudaMemoryManager->cudaAllocWallBC(lev);
+	  cudaMemoryManager->cudaAllocNoSlipBC(lev);
 
       findQ(para, lev);
  	  para->getParD(lev)->noSlipBC.numberOfBCnodes = para->getParH(lev)->noSlipBC.numberOfBCnodes;
       printf("numberOfBCnodes= %d\n", para->getParH(lev)->noSlipBC.numberOfBCnodes);
 
-	  cudaMemoryManager->cudaCopyWallBC(lev);
+	  cudaMemoryManager->cudaCopyNoSlipBC(lev);
    }
 }
 
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
index 5267f83838dcf131c3a6ab41ec05467e174ccfe8..a2026863e126d8648def4ffff6bf8dd22a4a0ba5 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
+++ b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
@@ -275,8 +275,8 @@ void CudaMemoryManager::cudaFreeOutflowBC(int lev)
 	checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->outflowBC.kN     ));
 	checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->outflowBC.RhoBC  ));
 }
-//Wall
-void CudaMemoryManager::cudaAllocWallBC(int lev)
+//No-Slip
+void CudaMemoryManager::cudaAllocNoSlipBC(int lev)
 {
 	unsigned int mem_size_Q_k      = sizeof(int)*parameter->getParH(lev)->noSlipBC.numberOfBCnodes;
 	unsigned int mem_size_Q_q      = sizeof(real)*parameter->getParH(lev)->noSlipBC.numberOfBCnodes;
@@ -297,7 +297,7 @@ void CudaMemoryManager::cudaAllocWallBC(int lev)
 	double tmp = (double)mem_size_Q_k + (double)parameter->getD3Qxx()*(double)mem_size_Q_q;
 	setMemsizeGPU(tmp, false);
 }
-void CudaMemoryManager::cudaCopyWallBC(int lev)
+void CudaMemoryManager::cudaCopyNoSlipBC(int lev)
 {
 	unsigned int mem_size_Q_k = sizeof(int)*parameter->getParH(lev)->noSlipBC.numberOfBCnodes;
 	unsigned int mem_size_Q_q = sizeof(real)*parameter->getParH(lev)->noSlipBC.numberOfBCnodes;
@@ -305,7 +305,7 @@ void CudaMemoryManager::cudaCopyWallBC(int lev)
 	checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->noSlipBC.q27[0], parameter->getParH(lev)->noSlipBC.q27[0], parameter->getD3Qxx()* mem_size_Q_q,       cudaMemcpyHostToDevice));
 	checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->noSlipBC.k,      parameter->getParH(lev)->noSlipBC.k,                  mem_size_Q_k,       cudaMemcpyHostToDevice));
 }
-void CudaMemoryManager::cudaFreeWallBC(int lev)
+void CudaMemoryManager::cudaFreeNoSlipBC(int lev)
 {
 	checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->noSlipBC.q27[0]));
 	checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->noSlipBC.k));
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h
index 290a05671e5f2d60b08d634ac2b919cd8947cb96..8df1b412c760ef3298f4426c206777d5c1db1b54 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h
+++ b/src/gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h
@@ -24,10 +24,10 @@ class Probe;
 class VIRTUALFLUIDS_GPU_EXPORT CudaMemoryManager
 {
 public:
-	CudaMemoryManager(std::shared_ptr<Parameter> parameter);
+    CudaMemoryManager(std::shared_ptr<Parameter> parameter);
 
-	void setMemsizeGPU(double admem, bool reset);
-	double getMemsizeGPU();
+    void setMemsizeGPU(double admem, bool reset);
+    double getMemsizeGPU();
 
     void cudaAllocFull(int lev);
     void cudaFreeFull(int lev);
@@ -35,109 +35,109 @@ public:
     void cudaCopyPrint(int lev);
     void cudaCopyMedianPrint(int lev);
 
-	void cudaAllocCoord(int lev);
-	void cudaCopyCoord(int lev);
-	void cudaFreeCoord(int lev);
+    void cudaAllocCoord(int lev);
+    void cudaCopyCoord(int lev);
+    void cudaFreeCoord(int lev);
 
-	void cudaAllocBodyForce(int lev);
+    void cudaAllocBodyForce(int lev);
     void cudaCopyBodyForce(int lev);
     void cudaFreeBodyForce(int lev);
 
     void cudaCopyDataToHost(int lev);
 
-	void cudaAllocSP(int lev);
-	void cudaCopySP(int lev);
-	void cudaFreeSP(int lev);
+    void cudaAllocSP(int lev);
+    void cudaCopySP(int lev);
+    void cudaFreeSP(int lev);
 
     void cudaAllocF3SP(int lev);
 
-	void cudaAllocNeighborWSB(int lev);
-	void cudaCopyNeighborWSB(int lev);
-	void cudaFreeNeighborWSB(int lev);
+    void cudaAllocNeighborWSB(int lev);
+    void cudaCopyNeighborWSB(int lev);
+    void cudaFreeNeighborWSB(int lev);
 
-	void cudaAllocVeloBC(int lev);
-	void cudaCopyVeloBC(int lev);
-	void cudaFreeVeloBC(int lev);
+    void cudaAllocVeloBC(int lev);
+    void cudaCopyVeloBC(int lev);
+    void cudaFreeVeloBC(int lev);
 
-	void cudaAllocOutflowBC(int lev);
-	void cudaCopyOutflowBC(int lev);
-	void cudaFreeOutflowBC(int lev);
+    void cudaAllocOutflowBC(int lev);
+    void cudaCopyOutflowBC(int lev);
+    void cudaFreeOutflowBC(int lev);
 
-	void cudaAllocWallBC(int lev);
-	void cudaCopyWallBC(int lev);
-	void cudaFreeWallBC(int lev);
+    void cudaAllocNoSlipBC(int lev);
+    void cudaCopyNoSlipBC(int lev);
+    void cudaFreeNoSlipBC(int lev);
 
-	void cudaAllocGeomBC(int lev);
-	void cudaCopyGeomBC(int lev);
-	void cudaFreeGeomBC(int lev);
+    void cudaAllocGeomBC(int lev);
+    void cudaCopyGeomBC(int lev);
+    void cudaFreeGeomBC(int lev);
 
-	void cudaAllocPress(int lev);
-	void cudaCopyPress(int lev);
-	void cudaFreePress(int lev);
+    void cudaAllocPress(int lev);
+    void cudaCopyPress(int lev);
+    void cudaFreePress(int lev);
 
-	void cudaAllocForcing();
-	void cudaCopyForcingToDevice();
-	void cudaCopyForcingToHost();
-	void cudaFreeForcing();
+    void cudaAllocForcing();
+    void cudaCopyForcingToDevice();
+    void cudaCopyForcingToHost();
+    void cudaFreeForcing();
 
     void cudaAllocLevelForcing(int level);
-	void cudaCopyLevelForcingToDevice(int level);
-	void cudaFreeLevelForcing(int level);
+    void cudaCopyLevelForcingToDevice(int level);
+    void cudaFreeLevelForcing(int level);
 
-	void cudaAllocQuadricLimiters();
-	void cudaCopyQuadricLimitersToDevice();
-	void cudaFreeQuadricLimiters();
+    void cudaAllocQuadricLimiters();
+    void cudaCopyQuadricLimitersToDevice();
+    void cudaFreeQuadricLimiters();
 
-	//////////////////////////////////////////////////////////////////////////
-	//3D domain decomposition
-	void cudaAllocProcessNeighborX(int lev, unsigned int processNeighbor);
+    //////////////////////////////////////////////////////////////////////////
+    //3D domain decomposition
+    void cudaAllocProcessNeighborX(int lev, unsigned int processNeighbor);
     void cudaCopyProcessNeighborXFsHD(int lev, unsigned int processNeighbor, const unsigned int &memsizeFsRecv,
                                       int streamIndex);
     void cudaCopyProcessNeighborXFsDH(int lev, unsigned int processNeighbor, const unsigned int &memsizeFsSend,
                                       int streamIndex);
-	void cudaCopyProcessNeighborXIndex(int lev, unsigned int processNeighbor);
-	void cudaFreeProcessNeighborX(int lev, unsigned int processNeighbor);
-	//
-	void cudaAllocProcessNeighborY(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborXIndex(int lev, unsigned int processNeighbor);
+    void cudaFreeProcessNeighborX(int lev, unsigned int processNeighbor);
+    //
+    void cudaAllocProcessNeighborY(int lev, unsigned int processNeighbor);
     void cudaCopyProcessNeighborYFsHD(int lev, unsigned int processNeighbor, const unsigned int &memsizeFsRecv,
                                       int streamIndex);
     void cudaCopyProcessNeighborYFsDH(int lev, unsigned int processNeighbor, const unsigned int &memsizeFsSend,
                                       int streamIndex);
-	void cudaCopyProcessNeighborYIndex(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborYIndex(int lev, unsigned int processNeighbor);
     void cudaFreeProcessNeighborY(int lev, unsigned int processNeighbor);
-	//
-	void cudaAllocProcessNeighborZ(int lev, unsigned int processNeighbor);
+    //
+    void cudaAllocProcessNeighborZ(int lev, unsigned int processNeighbor);
     void cudaCopyProcessNeighborZFsHD(int lev, unsigned int processNeighbor, const unsigned int &memsizeFsRecv,
                                       int streamIndex);
     void cudaCopyProcessNeighborZFsDH(int lev, unsigned int processNeighbor, const unsigned int &memsizeFsSend,
                                       int streamIndex);
-	void cudaCopyProcessNeighborZIndex(int lev, unsigned int processNeighbor);
-	void cudaFreeProcessNeighborZ(int lev, unsigned int processNeighbor);
-
-	//////////////////////////////////////////////////////////////////////////
-
-	//////////////////////////////////////////////////////////////////////////
-	//3D domain decomposition F3
-	void cudaAllocProcessNeighborF3X(int lev, unsigned int processNeighbor);
-	void cudaCopyProcessNeighborF3XFsHD(int lev, unsigned int processNeighbor);
-	void cudaCopyProcessNeighborF3XFsDH(int lev, unsigned int processNeighbor);
-	void cudaCopyProcessNeighborF3XIndex(int lev, unsigned int processNeighbor);
-	void cudaFreeProcessNeighborF3X(int lev, unsigned int processNeighbor);
-	//
-	void cudaAllocProcessNeighborF3Y(int lev, unsigned int processNeighbor);
-	void cudaCopyProcessNeighborF3YFsHD(int lev, unsigned int processNeighbor);
-	void cudaCopyProcessNeighborF3YFsDH(int lev, unsigned int processNeighbor);
-	void cudaCopyProcessNeighborF3YIndex(int lev, unsigned int processNeighbor);
-	void cudaFreeProcessNeighborF3Y(int lev, unsigned int processNeighbor);
-	//
-	void cudaAllocProcessNeighborF3Z(int lev, unsigned int processNeighbor);
-	void cudaCopyProcessNeighborF3ZFsHD(int lev, unsigned int processNeighbor);
-	void cudaCopyProcessNeighborF3ZFsDH(int lev, unsigned int processNeighbor);
-	void cudaCopyProcessNeighborF3ZIndex(int lev, unsigned int processNeighbor);
-	void cudaFreeProcessNeighborF3Z(int lev, unsigned int processNeighbor);
-	//////////////////////////////////////////////////////////////////////////
-
-	void cudaAllocTurbulentViscosity(int lev);
+    void cudaCopyProcessNeighborZIndex(int lev, unsigned int processNeighbor);
+    void cudaFreeProcessNeighborZ(int lev, unsigned int processNeighbor);
+
+    //////////////////////////////////////////////////////////////////////////
+
+    //////////////////////////////////////////////////////////////////////////
+    //3D domain decomposition F3
+    void cudaAllocProcessNeighborF3X(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborF3XFsHD(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborF3XFsDH(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborF3XIndex(int lev, unsigned int processNeighbor);
+    void cudaFreeProcessNeighborF3X(int lev, unsigned int processNeighbor);
+    //
+    void cudaAllocProcessNeighborF3Y(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborF3YFsHD(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborF3YFsDH(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborF3YIndex(int lev, unsigned int processNeighbor);
+    void cudaFreeProcessNeighborF3Y(int lev, unsigned int processNeighbor);
+    //
+    void cudaAllocProcessNeighborF3Z(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborF3ZFsHD(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborF3ZFsDH(int lev, unsigned int processNeighbor);
+    void cudaCopyProcessNeighborF3ZIndex(int lev, unsigned int processNeighbor);
+    void cudaFreeProcessNeighborF3Z(int lev, unsigned int processNeighbor);
+    //////////////////////////////////////////////////////////////////////////
+
+    void cudaAllocTurbulentViscosity(int lev);
     void cudaCopyTurbulentViscosityHD(int lev);
     void cudaCopyTurbulentViscosityDH(int lev);
     void cudaFreeTurbulentViscosity(int lev);
@@ -403,7 +403,7 @@ public:
 
 private:
     std::shared_ptr<Parameter> parameter;
-    double memsizeGPU = 0.;
+    double memsizeGPU = 0.0;
 
 };
 #endif
diff --git a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
index 7683d500b4360ed1660000e19f2c64f6e45ab4c2..316a49c420eaf76dcee49e23591bf166d1859f2e 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
@@ -3370,7 +3370,7 @@ void QSlipDevCompTurbulentViscosity27(LBMSimulationParameter* parameterDevice, Q
          parameterDevice->isEvenTimestep);
    getLastCudaError("QSlipDeviceComp27TurbViscosity execution failed");
 }
-
+//////////////////////////////////////////////////////////////////////////
 void QSlipDevComp27(LBMSimulationParameter* parameterDevice, QforBoundaryConditions* boundaryCondition)
 {
    dim3 grid = vf::cuda::getCudaGrid( parameterDevice->numberofthreads, boundaryCondition->numberOfBCnodes);
diff --git a/src/gpu/VirtualFluids_GPU/GPU/NoSlipBCs27.cu b/src/gpu/VirtualFluids_GPU/GPU/NoSlipBCs27.cu
index 4d720190aa21eb88abff99c3341cb7ad8e8a721d..31e83f56a6dd123d2d303374d6cf93adac9a32cd 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/NoSlipBCs27.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/NoSlipBCs27.cu
@@ -1641,7 +1641,7 @@ __global__ void QDeviceComp27(
 										 bool isEvenTimestep)
 {
    //////////////////////////////////////////////////////////////////////////
-   //! The velocity boundary condition is executed in the following steps
+   //! The no-slip boundary condition is executed in the following steps
    //!
    ////////////////////////////////////////////////////////////////////////////////
    //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
diff --git a/src/gpu/VirtualFluids_GPU/GPU/SlipBCs27.cu b/src/gpu/VirtualFluids_GPU/GPU/SlipBCs27.cu
index 1644e7c28ae700c8d1507fcec1f029ca463eeed0..717a41465ddedd5c170fddc273f7916b24f8cf49 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/SlipBCs27.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/SlipBCs27.cu
@@ -1,7 +1,8 @@
 /* Device code */
 #include "LBM/LB.h" 
 #include "lbm/constants/D3Q27.h"
-#include <lbm/constants/NumericConstants.h>
+#include "lbm/constants/NumericConstants.h"
+#include "KernelUtilities.h"
 
 using namespace vf::lbm::constant;
 using namespace vf::lbm::dir;
@@ -657,143 +658,68 @@ __global__ void QSlipDevice27(real* DD,
 
 
 //////////////////////////////////////////////////////////////////////////////
-__global__ void QSlipDeviceComp27(real* DD, 
-											 int* k_Q, 
-											 real* QQ,
-											 unsigned int numberOfBCnodes,
-											 real om1, 
-											 unsigned int* neighborX,
-											 unsigned int* neighborY,
-											 unsigned int* neighborZ,
-											 unsigned int size_Mat, 
-											 bool isEvenTimestep)
+__global__ void QSlipDeviceComp27(
+                                    real* distributions, 
+                                    int* subgridDistanceIndices, 
+                                    real* subgridDistances,
+                                    unsigned int numberOfBCnodes,
+                                    real omega, 
+                                    unsigned int* neighborX,
+                                    unsigned int* neighborY,
+                                    unsigned int* neighborZ,
+                                    unsigned int numberOfLBnodes, 
+                                    bool isEvenTimestep)
 {
-   Distributions27 D;
-   if (isEvenTimestep==true)
-   {
-      D.f[E   ] = &DD[E   *size_Mat];
-      D.f[W   ] = &DD[W   *size_Mat];
-      D.f[N   ] = &DD[N   *size_Mat];
-      D.f[S   ] = &DD[S   *size_Mat];
-      D.f[T   ] = &DD[T   *size_Mat];
-      D.f[B   ] = &DD[B   *size_Mat];
-      D.f[NE  ] = &DD[NE  *size_Mat];
-      D.f[SW  ] = &DD[SW  *size_Mat];
-      D.f[SE  ] = &DD[SE  *size_Mat];
-      D.f[NW  ] = &DD[NW  *size_Mat];
-      D.f[TE  ] = &DD[TE  *size_Mat];
-      D.f[BW  ] = &DD[BW  *size_Mat];
-      D.f[BE  ] = &DD[BE  *size_Mat];
-      D.f[TW  ] = &DD[TW  *size_Mat];
-      D.f[TN  ] = &DD[TN  *size_Mat];
-      D.f[BS  ] = &DD[BS  *size_Mat];
-      D.f[BN  ] = &DD[BN  *size_Mat];
-      D.f[TS  ] = &DD[TS  *size_Mat];
-      D.f[REST] = &DD[REST*size_Mat];
-      D.f[TNE ] = &DD[TNE *size_Mat];
-      D.f[TSW ] = &DD[TSW *size_Mat];
-      D.f[TSE ] = &DD[TSE *size_Mat];
-      D.f[TNW ] = &DD[TNW *size_Mat];
-      D.f[BNE ] = &DD[BNE *size_Mat];
-      D.f[BSW ] = &DD[BSW *size_Mat];
-      D.f[BSE ] = &DD[BSE *size_Mat];
-      D.f[BNW ] = &DD[BNW *size_Mat];
-   } 
-   else
-   {
-      D.f[W   ] = &DD[E   *size_Mat];
-      D.f[E   ] = &DD[W   *size_Mat];
-      D.f[S   ] = &DD[N   *size_Mat];
-      D.f[N   ] = &DD[S   *size_Mat];
-      D.f[B   ] = &DD[T   *size_Mat];
-      D.f[T   ] = &DD[B   *size_Mat];
-      D.f[SW  ] = &DD[NE  *size_Mat];
-      D.f[NE  ] = &DD[SW  *size_Mat];
-      D.f[NW  ] = &DD[SE  *size_Mat];
-      D.f[SE  ] = &DD[NW  *size_Mat];
-      D.f[BW  ] = &DD[TE  *size_Mat];
-      D.f[TE  ] = &DD[BW  *size_Mat];
-      D.f[TW  ] = &DD[BE  *size_Mat];
-      D.f[BE  ] = &DD[TW  *size_Mat];
-      D.f[BS  ] = &DD[TN  *size_Mat];
-      D.f[TN  ] = &DD[BS  *size_Mat];
-      D.f[TS  ] = &DD[BN  *size_Mat];
-      D.f[BN  ] = &DD[TS  *size_Mat];
-      D.f[REST] = &DD[REST*size_Mat];
-      D.f[TNE ] = &DD[BSW *size_Mat];
-      D.f[TSW ] = &DD[BNE *size_Mat];
-      D.f[TSE ] = &DD[BNW *size_Mat];
-      D.f[TNW ] = &DD[BSE *size_Mat];
-      D.f[BNE ] = &DD[TSW *size_Mat];
-      D.f[BSW ] = &DD[TNE *size_Mat];
-      D.f[BSE ] = &DD[TNW *size_Mat];
-      D.f[BNW ] = &DD[TSE *size_Mat];
-   }
+   //! The slip boundary condition is executed in the following steps
+   //!
    ////////////////////////////////////////////////////////////////////////////////
-   const unsigned  x = threadIdx.x;  // Globaler x-Index 
-   const unsigned  y = blockIdx.x;   // Globaler y-Index 
-   const unsigned  z = blockIdx.y;   // Globaler z-Index 
+   //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
+   //!
+   const unsigned  x = threadIdx.x;  // global x-index 
+   const unsigned  y = blockIdx.x;   // global y-index 
+   const unsigned  z = blockIdx.y;   // global z-index 
 
    const unsigned nx = blockDim.x;
    const unsigned ny = gridDim.x;
 
    const unsigned k = nx*(ny*z + y) + x;
-   //////////////////////////////////////////////////////////////////////////
 
    if(k < numberOfBCnodes)
    {
+      //////////////////////////////////////////////////////////////////////////
+      //! - Read distributions: style of reading and writing the distributions from/to stored arrays dependent on timestep is based on the esoteric twist algorithm \ref
+      //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
+      //!
+      Distributions27 dist;
+      getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
+
       ////////////////////////////////////////////////////////////////////////////////
-      real *q_dirE,   *q_dirW,   *q_dirN,   *q_dirS,   *q_dirT,   *q_dirB, 
-            *q_dirNE,  *q_dirSW,  *q_dirSE,  *q_dirNW,  *q_dirTE,  *q_dirBW,
-            *q_dirBE,  *q_dirTW,  *q_dirTN,  *q_dirBS,  *q_dirBN,  *q_dirTS,
-            *q_dirTNE, *q_dirTSW, *q_dirTSE, *q_dirTNW, *q_dirBNE, *q_dirBSW,
-            *q_dirBSE, *q_dirBNW; 
-      q_dirE   = &QQ[E   * numberOfBCnodes];
-      q_dirW   = &QQ[W   * numberOfBCnodes];
-      q_dirN   = &QQ[N   * numberOfBCnodes];
-      q_dirS   = &QQ[S   * numberOfBCnodes];
-      q_dirT   = &QQ[T   * numberOfBCnodes];
-      q_dirB   = &QQ[B   * numberOfBCnodes];
-      q_dirNE  = &QQ[NE  * numberOfBCnodes];
-      q_dirSW  = &QQ[SW  * numberOfBCnodes];
-      q_dirSE  = &QQ[SE  * numberOfBCnodes];
-      q_dirNW  = &QQ[NW  * numberOfBCnodes];
-      q_dirTE  = &QQ[TE  * numberOfBCnodes];
-      q_dirBW  = &QQ[BW  * numberOfBCnodes];
-      q_dirBE  = &QQ[BE  * numberOfBCnodes];
-      q_dirTW  = &QQ[TW  * numberOfBCnodes];
-      q_dirTN  = &QQ[TN  * numberOfBCnodes];
-      q_dirBS  = &QQ[BS  * numberOfBCnodes];
-      q_dirBN  = &QQ[BN  * numberOfBCnodes];
-      q_dirTS  = &QQ[TS  * numberOfBCnodes];
-      q_dirTNE = &QQ[TNE * numberOfBCnodes];
-      q_dirTSW = &QQ[TSW * numberOfBCnodes];
-      q_dirTSE = &QQ[TSE * numberOfBCnodes];
-      q_dirTNW = &QQ[TNW * numberOfBCnodes];
-      q_dirBNE = &QQ[BNE * numberOfBCnodes];
-      q_dirBSW = &QQ[BSW * numberOfBCnodes];
-      q_dirBSE = &QQ[BSE * numberOfBCnodes];
-      q_dirBNW = &QQ[BNW * numberOfBCnodes];
+      //! - Set local subgrid distances (q's)
+      //!
+      SubgridDistances27 subgridD;
+      getPointersToSubgridDistances(subgridD, subgridDistances, numberOfBCnodes);
+      
       ////////////////////////////////////////////////////////////////////////////////
-      //index
-      unsigned int KQK  = k_Q[k];
-      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];
+      //! - Set neighbor indices (necessary for indirect addressing)
+      //!
+      unsigned int indexOfBCnode  = subgridDistanceIndices[k];
+      unsigned int kzero= indexOfBCnode;
+      unsigned int ke   = indexOfBCnode;
+      unsigned int kw   = neighborX[indexOfBCnode];
+      unsigned int kn   = indexOfBCnode;
+      unsigned int ks   = neighborY[indexOfBCnode];
+      unsigned int kt   = indexOfBCnode;
+      unsigned int kb   = neighborZ[indexOfBCnode];
       unsigned int ksw  = neighborY[kw];
-      unsigned int kne  = KQK;
+      unsigned int kne  = indexOfBCnode;
       unsigned int kse  = ks;
       unsigned int knw  = kw;
       unsigned int kbw  = neighborZ[kw];
-      unsigned int kte  = KQK;
+      unsigned int kte  = indexOfBCnode;
       unsigned int kbe  = kb;
       unsigned int ktw  = kw;
       unsigned int kbs  = neighborZ[ks];
-      unsigned int ktn  = KQK;
+      unsigned int ktn  = indexOfBCnode;
       unsigned int kbn  = kb;
       unsigned int kts  = ks;
       unsigned int ktse = ks;
@@ -802,521 +728,442 @@ __global__ void QSlipDeviceComp27(real* DD,
       unsigned int kbse = kbs;
       unsigned int ktsw = ksw;
       unsigned int kbne = kb;
-      unsigned int ktne = KQK;
+      unsigned int ktne = indexOfBCnode;
       unsigned int kbsw = neighborZ[ksw];
       
       ////////////////////////////////////////////////////////////////////////////////
-      real f_W    = (D.f[E   ])[ke   ];
-      real f_E    = (D.f[W   ])[kw   ];
-      real f_S    = (D.f[N   ])[kn   ];
-      real f_N    = (D.f[S   ])[ks   ];
-      real f_B    = (D.f[T   ])[kt   ];
-      real f_T    = (D.f[B   ])[kb   ];
-      real f_SW   = (D.f[NE  ])[kne  ];
-      real f_NE   = (D.f[SW  ])[ksw  ];
-      real f_NW   = (D.f[SE  ])[kse  ];
-      real f_SE   = (D.f[NW  ])[knw  ];
-      real f_BW   = (D.f[TE  ])[kte  ];
-      real f_TE   = (D.f[BW  ])[kbw  ];
-      real f_TW   = (D.f[BE  ])[kbe  ];
-      real f_BE   = (D.f[TW  ])[ktw  ];
-      real f_BS   = (D.f[TN  ])[ktn  ];
-      real f_TN   = (D.f[BS  ])[kbs  ];
-      real f_TS   = (D.f[BN  ])[kbn  ];
-      real f_BN   = (D.f[TS  ])[kts  ];
-      real f_BSW  = (D.f[TNE ])[ktne ];
-      real f_BNE  = (D.f[TSW ])[ktsw ];
-      real f_BNW  = (D.f[TSE ])[ktse ];
-      real f_BSE  = (D.f[TNW ])[ktnw ];
-      real f_TSW  = (D.f[BNE ])[kbne ];
-      real f_TNE  = (D.f[BSW ])[kbsw ];
-      real f_TNW  = (D.f[BSE ])[kbse ];
-      real f_TSE  = (D.f[BNW ])[kbnw ];
+      //! - Set local distributions
+      //!
+      real f_W    = (dist.f[E   ])[ke   ];
+      real f_E    = (dist.f[W   ])[kw   ];
+      real f_S    = (dist.f[N   ])[kn   ];
+      real f_N    = (dist.f[S   ])[ks   ];
+      real f_B    = (dist.f[T   ])[kt   ];
+      real f_T    = (dist.f[B   ])[kb   ];
+      real f_SW   = (dist.f[NE  ])[kne  ];
+      real f_NE   = (dist.f[SW  ])[ksw  ];
+      real f_NW   = (dist.f[SE  ])[kse  ];
+      real f_SE   = (dist.f[NW  ])[knw  ];
+      real f_BW   = (dist.f[TE  ])[kte  ];
+      real f_TE   = (dist.f[BW  ])[kbw  ];
+      real f_TW   = (dist.f[BE  ])[kbe  ];
+      real f_BE   = (dist.f[TW  ])[ktw  ];
+      real f_BS   = (dist.f[TN  ])[ktn  ];
+      real f_TN   = (dist.f[BS  ])[kbs  ];
+      real f_TS   = (dist.f[BN  ])[kbn  ];
+      real f_BN   = (dist.f[TS  ])[kts  ];
+      real f_BSW  = (dist.f[TNE ])[ktne ];
+      real f_BNE  = (dist.f[TSW ])[ktsw ];
+      real f_BNW  = (dist.f[TSE ])[ktse ];
+      real f_BSE  = (dist.f[TNW ])[ktnw ];
+      real f_TSW  = (dist.f[BNE ])[kbne ];
+      real f_TNE  = (dist.f[BSW ])[kbsw ];
+      real f_TNW  = (dist.f[BSE ])[kbse ];
+      real f_TSE  = (dist.f[BNW ])[kbnw ];
+
       ////////////////////////////////////////////////////////////////////////////////
-      real vx1, vx2, vx3, drho, feq, q;
-      drho   =  f_TSE + f_TNW + f_TNE + f_TSW + f_BSE + f_BNW + f_BNE + f_BSW +
-                f_BN + f_TS + f_TN + f_BS + f_BE + f_TW + f_TE + f_BW + f_SE + f_NW + f_NE + f_SW + 
-                f_T + f_B + f_N + f_S + f_E + f_W + ((D.f[REST])[kzero]); 
+      //! - Calculate macroscopic quantities
+      //!
+      real drho = f_TSE + f_TNW + f_TNE + f_TSW + f_BSE + f_BNW + f_BNE + f_BSW +
+                  f_BN + f_TS + f_TN + f_BS + f_BE + f_TW + f_TE + f_BW + f_SE + f_NW + f_NE + f_SW + 
+                  f_T + f_B + f_N + f_S + f_E + f_W + ((dist.f[REST])[kzero]); 
 
-      vx1    =  (((f_TSE - f_BNW) - (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
-                ((f_BE - f_TW)   + (f_TE - f_BW))   + ((f_SE - f_NW)   + (f_NE - f_SW)) +
-                (f_E - f_W)) / (c1o1 + drho); 
-         
+      real vx1  = (((f_TSE - f_BNW) - (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
+                   ((f_BE - f_TW)   + (f_TE - f_BW))   + ((f_SE - f_NW)   + (f_NE - f_SW)) +
+                   (f_E - f_W)) / (c1o1 + drho);
 
-      vx2    =   ((-(f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
-                 ((f_BN - f_TS)   + (f_TN - f_BS))    + (-(f_SE - f_NW)  + (f_NE - f_SW)) +
-                 (f_N - f_S)) / (c1o1 + drho); 
+      real vx2  = ((-(f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
+                   ((f_BN - f_TS)   + (f_TN - f_BS))    + (-(f_SE - f_NW)  + (f_NE - f_SW)) +
+                   (f_N - f_S)) / (c1o1 + drho);
 
-      vx3    =   (((f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) + (f_TSW - f_BNE)) +
-                 (-(f_BN - f_TS)  + (f_TN - f_BS))   + ((f_TE - f_BW)   - (f_BE - f_TW)) +
-                 (f_T - f_B)) / (c1o1 + drho); 
+      real vx3  = (((f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) + (f_TSW - f_BNE)) +
+                   (-(f_BN - f_TS)  + (f_TN - f_BS))   + ((f_TE - f_BW)   - (f_BE - f_TW)) +
+                   (f_T - f_B)) / (c1o1 + drho);
 
-      real cu_sq=c3o2*(vx1*vx1+vx2*vx2+vx3*vx3) * (c1o1 + drho);
+      real cu_sq = c3o2 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3) * (c1o1 + drho);
 
-      //////////////////////////////////////////////////////////////////////////
-      if (isEvenTimestep==false)
-      {
-         D.f[E   ] = &DD[E   *size_Mat];
-         D.f[W   ] = &DD[W   *size_Mat];
-         D.f[N   ] = &DD[N   *size_Mat];
-         D.f[S   ] = &DD[S   *size_Mat];
-         D.f[T   ] = &DD[T   *size_Mat];
-         D.f[B   ] = &DD[B   *size_Mat];
-         D.f[NE  ] = &DD[NE  *size_Mat];
-         D.f[SW  ] = &DD[SW  *size_Mat];
-         D.f[SE  ] = &DD[SE  *size_Mat];
-         D.f[NW  ] = &DD[NW  *size_Mat];
-         D.f[TE  ] = &DD[TE  *size_Mat];
-         D.f[BW  ] = &DD[BW  *size_Mat];
-         D.f[BE  ] = &DD[BE  *size_Mat];
-         D.f[TW  ] = &DD[TW  *size_Mat];
-         D.f[TN  ] = &DD[TN  *size_Mat];
-         D.f[BS  ] = &DD[BS  *size_Mat];
-         D.f[BN  ] = &DD[BN  *size_Mat];
-         D.f[TS  ] = &DD[TS  *size_Mat];
-         D.f[REST] = &DD[REST*size_Mat];
-         D.f[TNE ] = &DD[TNE *size_Mat];
-         D.f[TSW ] = &DD[TSW *size_Mat];
-         D.f[TSE ] = &DD[TSE *size_Mat];
-         D.f[TNW ] = &DD[TNW *size_Mat];
-         D.f[BNE ] = &DD[BNE *size_Mat];
-         D.f[BSW ] = &DD[BSW *size_Mat];
-         D.f[BSE ] = &DD[BSE *size_Mat];
-         D.f[BNW ] = &DD[BNW *size_Mat];
-      } 
-      else
-      {
-         D.f[W   ] = &DD[E   *size_Mat];
-         D.f[E   ] = &DD[W   *size_Mat];
-         D.f[S   ] = &DD[N   *size_Mat];
-         D.f[N   ] = &DD[S   *size_Mat];
-         D.f[B   ] = &DD[T   *size_Mat];
-         D.f[T   ] = &DD[B   *size_Mat];
-         D.f[SW  ] = &DD[NE  *size_Mat];
-         D.f[NE  ] = &DD[SW  *size_Mat];
-         D.f[NW  ] = &DD[SE  *size_Mat];
-         D.f[SE  ] = &DD[NW  *size_Mat];
-         D.f[BW  ] = &DD[TE  *size_Mat];
-         D.f[TE  ] = &DD[BW  *size_Mat];
-         D.f[TW  ] = &DD[BE  *size_Mat];
-         D.f[BE  ] = &DD[TW  *size_Mat];
-         D.f[BS  ] = &DD[TN  *size_Mat];
-         D.f[TN  ] = &DD[BS  *size_Mat];
-         D.f[TS  ] = &DD[BN  *size_Mat];
-         D.f[BN  ] = &DD[TS  *size_Mat];
-         D.f[REST] = &DD[REST*size_Mat];
-         D.f[TNE ] = &DD[BSW *size_Mat];
-         D.f[TSW ] = &DD[BNE *size_Mat];
-         D.f[TSE ] = &DD[BNW *size_Mat];
-         D.f[TNW ] = &DD[BSE *size_Mat];
-         D.f[BNE ] = &DD[TSW *size_Mat];
-         D.f[BSW ] = &DD[TNE *size_Mat];
-         D.f[BSE ] = &DD[TNW *size_Mat];
-         D.f[BNW ] = &DD[TSE *size_Mat];
-      }
-      ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-      //Test
-      //(D.f[REST])[k]=c1o10;
-      ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-	  real fac = c1o1;//c99o100;
-	  real VeloX = fac*vx1;
-	  real VeloY = fac*vx2;
-	  real VeloZ = fac*vx3;
-	  bool x = false;
-	  bool y = false;
-	  bool z = false;
+      ////////////////////////////////////////////////////////////////////////////////
+      //! - change the pointer to write the results in the correct array
+      //!
+      getPointersToDistributions(dist, distributions, numberOfLBnodes, !isEvenTimestep);
 
-      q = q_dirE[k];
-      if (q>=c0o1 && q<=c1o1)
+      ////////////////////////////////////////////////////////////////////////////////
+      //! - Multiply the local velocities by the slipLength
+      //!
+      real slipLength = c1o1;
+      real VeloX = slipLength*vx1;
+      real VeloY = slipLength*vx2;
+      real VeloZ = slipLength*vx3;
+
+      ////////////////////////////////////////////////////////////////////////////////
+      //! - Update distributions with subgrid distance (q) between zero and one
+      //!
+      real feq, q, velocityLB, velocityBC;
+
+      bool x = false;
+      bool y = false;
+      bool z = false;
+
+      q = (subgridD.q[E])[k];
+      if (q>=c0o1 && q<=c1o1)  // only update distribution for q between zero and one
       {
-		 VeloX = c0o1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 x = true;
-         feq=c2o27* (drho/*+three*( vx1        )*/+c9o2*( vx1        )*( vx1        ) * (c1o1 + drho)-cu_sq); 
-         (D.f[W])[kw]=(c1o1-q)/(c1o1+q)*(f_E-f_W+(f_E+f_W-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_E+f_W)-c6o1*c2o27*( VeloX     ))/(c1o1+q) - c2o27 * drho;
-         //feq=c2over27* (drho+three*( vx1        )+c9over2*( vx1        )*( vx1        )-cu_sq); 
-         //(D.f[W])[kw]=(one-q)/(one+q)*(f_E-feq*om1)/(one-om1)+(q*(f_E+f_W)-six*c2over27*( VeloX     ))/(one+q);
-         //(D.f[W])[kw]=zero;
+         VeloX = c0o1;
+         x = true;
+
+         velocityLB = vx1;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
+         velocityBC = VeloX;
+         (dist.f[W])[kw] = getInterpolatedDistributionForVeloBC(q, f_E, f_W, feq, omega, velocityBC, c2o27);
       }
 
-      q = q_dirW[k];
+      q = (subgridD.q[W])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = c0o1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 x = true;
-         feq=c2o27* (drho/*+three*(-vx1        )*/+c9o2*(-vx1        )*(-vx1        ) * (c1o1 + drho)-cu_sq); 
-         (D.f[E])[ke]=(c1o1-q)/(c1o1+q)*(f_W-f_E+(f_W+f_E-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_W+f_E)-c6o1*c2o27*(-VeloX     ))/(c1o1+q) - c2o27 * drho;
-         //feq=c2over27* (drho+three*(-vx1        )+c9over2*(-vx1        )*(-vx1        )-cu_sq); 
-         //(D.f[E])[ke]=(one-q)/(one+q)*(f_W-feq*om1)/(one-om1)+(q*(f_W+f_E)-six*c2over27*(-VeloX     ))/(one+q);
-         //(D.f[E])[ke]=zero;
+         VeloX = c0o1;
+         x = true;
+
+         velocityLB = -vx1;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
+         velocityBC = -VeloX;
+         (dist.f[E])[ke] = getInterpolatedDistributionForVeloBC(q, f_W, f_E, feq, omega, velocityBC, c2o27);
       }
 
-      q = q_dirN[k];
-      if (q>=c0o1 && q<=c1o1)
-      {
-		 VeloX = fac*vx1;
-		 VeloY = c0o1;
-	     VeloZ = fac*vx3;
-		 y = true;
-         feq=c2o27* (drho/*+three*(    vx2     )*/+c9o2*(     vx2    )*(     vx2    ) * (c1o1 + drho)-cu_sq); 
-         (D.f[S])[ks]=(c1o1-q)/(c1o1+q)*(f_N-f_S+(f_N+f_S-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_N+f_S)-c6o1*c2o27*( VeloY     ))/(c1o1+q) - c2o27 * drho;
-         //feq=c2over27* (drho+three*(    vx2     )+c9over2*(     vx2    )*(     vx2    )-cu_sq); 
-         //(D.f[S])[ks]=(one-q)/(one+q)*(f_N-feq*om1)/(one-om1)+(q*(f_N+f_S)-six*c2over27*( VeloY     ))/(one+q);
-         //(D.f[S])[ks]=zero;
+      q = (subgridD.q[N])[k];
+      if (q>=c0o1 && q<=c1o1)
+      {
+         VeloY = c0o1;
+         y = true;
+
+         velocityLB = vx2;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
+         velocityBC = VeloY;
+         (dist.f[S])[ks] = getInterpolatedDistributionForVeloBC(q, f_N, f_S, feq, omega, velocityBC, c2o27);
       }
 
-      q = q_dirS[k];
+      q = (subgridD.q[S])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-		 VeloY = c0o1;
-	     VeloZ = fac*vx3;
-		 y = true;
-         feq=c2o27* (drho/*+three*(   -vx2     )*/+c9o2*(    -vx2    )*(    -vx2    ) * (c1o1 + drho)-cu_sq); 
-         (D.f[N])[kn]=(c1o1-q)/(c1o1+q)*(f_S-f_N+(f_S+f_N-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_S+f_N)-c6o1*c2o27*(-VeloY     ))/(c1o1+q) - c2o27 * drho;
-         //feq=c2over27* (drho+three*(   -vx2     )+c9over2*(    -vx2    )*(    -vx2    )-cu_sq); 
-         //(D.f[N])[kn]=(one-q)/(one+q)*(f_S-feq*om1)/(one-om1)+(q*(f_S+f_N)-six*c2over27*(-VeloY     ))/(one+q);
-         //(D.f[N])[kn]=zero;
+         VeloY = c0o1;
+         y = true;
+
+         velocityLB = -vx2;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
+         velocityBC = -VeloY;
+         (dist.f[N])[kn] = getInterpolatedDistributionForVeloBC(q, f_S, f_N, feq, omega, velocityBC, c2o27);
       }
 
-      q = q_dirT[k];
+      q = (subgridD.q[T])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-		 VeloZ = c0o1;
-		 z = true;
-         feq=c2o27* (drho/*+three*(         vx3)*/+c9o2*(         vx3)*(         vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[B])[kb]=(c1o1-q)/(c1o1+q)*(f_T-f_B+(f_T+f_B-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_T+f_B)-c6o1*c2o27*( VeloZ     ))/(c1o1+q) - c2o27 * drho;
-         //feq=c2over27* (drho+three*(         vx3)+c9over2*(         vx3)*(         vx3)-cu_sq); 
-         //(D.f[B])[kb]=(one-q)/(one+q)*(f_T-feq*om1)/(one-om1)+(q*(f_T+f_B)-six*c2over27*( VeloZ     ))/(one+q);
-         //(D.f[B])[kb]=one;
+         VeloZ = c0o1;
+         z = true;
+
+         velocityLB = vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
+         velocityBC = VeloZ;
+         (dist.f[B])[kb] = getInterpolatedDistributionForVeloBC(q, f_T, f_B, feq, omega, velocityBC, c2o27);
       }
 
-      q = q_dirB[k];
+      q = (subgridD.q[B])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-		 VeloZ = c0o1;
-		 z = true;
-         feq=c2o27* (drho/*+three*(        -vx3)*/+c9o2*(        -vx3)*(        -vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[T])[kt]=(c1o1-q)/(c1o1+q)*(f_B-f_T+(f_B+f_T-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_B+f_T)-c6o1*c2o27*(-VeloZ     ))/(c1o1+q) - c2o27 * drho;
-         //feq=c2over27* (drho+three*(        -vx3)+c9over2*(        -vx3)*(        -vx3)-cu_sq); 
-         //(D.f[T])[kt]=(one-q)/(one+q)*(f_B-feq*om1)/(one-om1)+(q*(f_B+f_T)-six*c2over27*(-VeloZ     ))/(one+q);
-         //(D.f[T])[kt]=zero;
+         VeloZ = c0o1;
+         z = true;
+
+         velocityLB = -vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c2o27);
+         velocityBC = -VeloZ;
+         (dist.f[T])[kt] = getInterpolatedDistributionForVeloBC(q, f_B, f_T, feq, omega, velocityBC, c2o27);
       }
 
-      q = q_dirNE[k];
+      q = (subgridD.q[NE])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-         feq=c1o54* (drho/*+three*( vx1+vx2    )*/+c9o2*( vx1+vx2    )*( vx1+vx2    ) * (c1o1 + drho)-cu_sq); 
-         (D.f[SW])[ksw]=(c1o1-q)/(c1o1+q)*(f_NE-f_SW+(f_NE+f_SW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_NE+f_SW)-c6o1*c1o54*(VeloX+VeloY))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*( vx1+vx2    )+c9over2*( vx1+vx2    )*( vx1+vx2    )-cu_sq); 
-         //(D.f[SW])[ksw]=(one-q)/(one+q)*(f_NE-feq*om1)/(one-om1)+(q*(f_NE+f_SW)-six*c1over54*(VeloX+VeloY))/(one+q);
-         //(D.f[SW])[ksw]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+
+         velocityLB = vx1 + vx2;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = VeloX + VeloY;
+         (dist.f[SW])[ksw] = getInterpolatedDistributionForVeloBC(q, f_NE, f_SW, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirSW[k];
+      q = (subgridD.q[SW])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-         feq=c1o54* (drho/*+three*(-vx1-vx2    )*/+c9o2*(-vx1-vx2    )*(-vx1-vx2    ) * (c1o1 + drho)-cu_sq); 
-         (D.f[NE])[kne]=(c1o1-q)/(c1o1+q)*(f_SW-f_NE+(f_SW+f_NE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_SW+f_NE)-c6o1*c1o54*(-VeloX-VeloY))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*(-vx1-vx2    )+c9over2*(-vx1-vx2    )*(-vx1-vx2    )-cu_sq); 
-         //(D.f[NE])[kne]=(one-q)/(one+q)*(f_SW-feq*om1)/(one-om1)+(q*(f_SW+f_NE)-six*c1over54*(-VeloX-VeloY))/(one+q);
-         //(D.f[NE])[kne]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+
+         velocityLB = -vx1 - vx2;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = -VeloX - VeloY;
+         (dist.f[NE])[kne] = getInterpolatedDistributionForVeloBC(q, f_SW, f_NE, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirSE[k];
+      q = (subgridD.q[SE])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-         feq=c1o54* (drho/*+three*( vx1-vx2    )*/+c9o2*( vx1-vx2    )*( vx1-vx2    ) * (c1o1 + drho)-cu_sq); 
-         (D.f[NW])[knw]=(c1o1-q)/(c1o1+q)*(f_SE-f_NW+(f_SE+f_NW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_SE+f_NW)-c6o1*c1o54*( VeloX-VeloY))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*( vx1-vx2    )+c9over2*( vx1-vx2    )*( vx1-vx2    )-cu_sq); 
-         //(D.f[NW])[knw]=(one-q)/(one+q)*(f_SE-feq*om1)/(one-om1)+(q*(f_SE+f_NW)-six*c1over54*( VeloX-VeloY))/(one+q);
-         //(D.f[NW])[knw]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+
+         velocityLB = vx1 - vx2;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = VeloX - VeloY;
+         (dist.f[NW])[knw] = getInterpolatedDistributionForVeloBC(q, f_SE, f_NW, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirNW[k];
+      q = (subgridD.q[NW])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-         feq=c1o54* (drho/*+three*(-vx1+vx2    )*/+c9o2*(-vx1+vx2    )*(-vx1+vx2    ) * (c1o1 + drho)-cu_sq); 
-         (D.f[SE])[kse]=(c1o1-q)/(c1o1+q)*(f_NW-f_SE+(f_NW+f_SE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_NW+f_SE)-c6o1*c1o54*(-VeloX+VeloY))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*(-vx1+vx2    )+c9over2*(-vx1+vx2    )*(-vx1+vx2    )-cu_sq); 
-         //(D.f[SE])[kse]=(one-q)/(one+q)*(f_NW-feq*om1)/(one-om1)+(q*(f_NW+f_SE)-six*c1over54*(-VeloX+VeloY))/(one+q);
-         //(D.f[SE])[kse]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+
+         velocityLB = -vx1 + vx2;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = -VeloX + VeloY;
+         (dist.f[SE])[kse] = getInterpolatedDistributionForVeloBC(q, f_NW, f_SE, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirTE[k];
+      q = (subgridD.q[TE])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (z == true) VeloZ = c0o1;
-      //  if (k==10000) printf("AFTER x: %u \t  y: %u \t z: %u \n  VeloX: %f \t VeloY: %f \t VeloZ: %f \n\n", x,y,z, VeloX,VeloY,VeloZ);
-         feq=c1o54* (drho/*+three*( vx1    +vx3)*/+c9o2*( vx1    +vx3)*( vx1    +vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[BW])[kbw]=(c1o1-q)/(c1o1+q)*(f_TE-f_BW+(f_TE+f_BW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TE+f_BW)-c6o1*c1o54*( VeloX+VeloZ))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*( vx1    +vx3)+c9over2*( vx1    +vx3)*( vx1    +vx3)-cu_sq); 
-         //(D.f[BW])[kbw]=(one-q)/(one+q)*(f_TE-feq*om1)/(one-om1)+(q*(f_TE+f_BW)-six*c1over54*( VeloX+VeloZ))/(one+q);
-         //(D.f[BW])[kbw]=zero;
+         VeloX = slipLength*vx1;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (z == true) VeloZ = c0o1;
+
+         velocityLB = vx1 + vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = VeloX + VeloZ;
+         (dist.f[BW])[kbw] = getInterpolatedDistributionForVeloBC(q, f_TE, f_BW, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirBW[k];
+      q = (subgridD.q[BW])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o54* (drho/*+three*(-vx1    -vx3)*/+c9o2*(-vx1    -vx3)*(-vx1    -vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[TE])[kte]=(c1o1-q)/(c1o1+q)*(f_BW-f_TE+(f_BW+f_TE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BW+f_TE)-c6o1*c1o54*(-VeloX-VeloZ))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*(-vx1    -vx3)+c9over2*(-vx1    -vx3)*(-vx1    -vx3)-cu_sq); 
-         //(D.f[TE])[kte]=(one-q)/(one+q)*(f_BW-feq*om1)/(one-om1)+(q*(f_BW+f_TE)-six*c1over54*(-VeloX-VeloZ))/(one+q);
-         //(D.f[TE])[kte]=zero;
+        VeloX = slipLength*vx1;
+        VeloZ = slipLength*vx3;
+        if (x == true) VeloX = c0o1;
+        if (z == true) VeloZ = c0o1;
+
+         velocityLB = -vx1 - vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = -VeloX - VeloZ;
+         (dist.f[TE])[kte] = getInterpolatedDistributionForVeloBC(q, f_BW, f_TE, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirBE[k];
+      q = (subgridD.q[BE])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o54* (drho/*+three*( vx1    -vx3)*/+c9o2*( vx1    -vx3)*( vx1    -vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[TW])[ktw]=(c1o1-q)/(c1o1+q)*(f_BE-f_TW+(f_BE+f_TW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BE+f_TW)-c6o1*c1o54*( VeloX-VeloZ))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*( vx1    -vx3)+c9over2*( vx1    -vx3)*( vx1    -vx3)-cu_sq); 
-         //(D.f[TW])[ktw]=(one-q)/(one+q)*(f_BE-feq*om1)/(one-om1)+(q*(f_BE+f_TW)-six*c1over54*( VeloX-VeloZ))/(one+q);
-         //(D.f[TW])[ktw]=zero;
+         VeloX = slipLength*vx1;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (z == true) VeloZ = c0o1;
+
+         velocityLB = vx1 - vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = VeloX - VeloZ;
+         (dist.f[TW])[ktw] = getInterpolatedDistributionForVeloBC(q, f_BE, f_TW, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirTW[k];
+      q = (subgridD.q[TW])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o54* (drho/*+three*(-vx1    +vx3)*/+c9o2*(-vx1    +vx3)*(-vx1    +vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[BE])[kbe]=(c1o1-q)/(c1o1+q)*(f_TW-f_BE+(f_TW+f_BE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TW+f_BE)-c6o1*c1o54*(-VeloX+VeloZ))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*(-vx1    +vx3)+c9over2*(-vx1    +vx3)*(-vx1    +vx3)-cu_sq); 
-         //(D.f[BE])[kbe]=(one-q)/(one+q)*(f_TW-feq*om1)/(one-om1)+(q*(f_TW+f_BE)-six*c1over54*(-VeloX+VeloZ))/(one+q);
-         //(D.f[BE])[kbe]=zero;
+         VeloX = slipLength*vx1;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (z == true) VeloZ = c0o1;
+
+         velocityLB = -vx1 + vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = -VeloX + VeloZ;
+         (dist.f[BE])[kbe] = getInterpolatedDistributionForVeloBC(q, f_TW, f_BE, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirTN[k];
+      q = (subgridD.q[TN])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o54* (drho/*+three*(     vx2+vx3)*/+c9o2*(     vx2+vx3)*(     vx2+vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[BS])[kbs]=(c1o1-q)/(c1o1+q)*(f_TN-f_BS+(f_TN+f_BS-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TN+f_BS)-c6o1*c1o54*( VeloY+VeloZ))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*(     vx2+vx3)+c9over2*(     vx2+vx3)*(     vx2+vx3)-cu_sq); 
-         //(D.f[BS])[kbs]=(one-q)/(one+q)*(f_TN-feq*om1)/(one-om1)+(q*(f_TN+f_BS)-six*c1over54*( VeloY+VeloZ))/(one+q);
-         //(D.f[BS])[kbs]=zero;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+
+         velocityLB = vx2 + vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = VeloY + VeloZ;
+         (dist.f[BS])[kbs] = getInterpolatedDistributionForVeloBC(q, f_TN, f_BS, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirBS[k];
+      q = (subgridD.q[BS])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o54* (drho/*+three*(    -vx2-vx3)*/+c9o2*(    -vx2-vx3)*(    -vx2-vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[TN])[ktn]=(c1o1-q)/(c1o1+q)*(f_BS-f_TN+(f_BS+f_TN-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BS+f_TN)-c6o1*c1o54*( -VeloY-VeloZ))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*(    -vx2-vx3)+c9over2*(    -vx2-vx3)*(    -vx2-vx3)-cu_sq); 
-         //(D.f[TN])[ktn]=(one-q)/(one+q)*(f_BS-feq*om1)/(one-om1)+(q*(f_BS+f_TN)-six*c1over54*( -VeloY-VeloZ))/(one+q);
-         //(D.f[TN])[ktn]=zero;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+
+         velocityLB = -vx2 - vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = -VeloY - VeloZ;
+         (dist.f[TN])[ktn] = getInterpolatedDistributionForVeloBC(q, f_BS, f_TN, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirBN[k];
+
+      q = (subgridD.q[BN])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o54* (drho/*+three*(     vx2-vx3)*/+c9o2*(     vx2-vx3)*(     vx2-vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[TS])[kts]=(c1o1-q)/(c1o1+q)*(f_BN-f_TS+(f_BN+f_TS-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BN+f_TS)-c6o1*c1o54*( VeloY-VeloZ))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*(     vx2-vx3)+c9over2*(     vx2-vx3)*(     vx2-vx3)-cu_sq); 
-         //(D.f[TS])[kts]=(one-q)/(one+q)*(f_BN-feq*om1)/(one-om1)+(q*(f_BN+f_TS)-six*c1over54*( VeloY-VeloZ))/(one+q);
-         //(D.f[TS])[kts]=zero;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+
+         velocityLB = vx2 - vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = VeloY - VeloZ;
+         (dist.f[TS])[kts] = getInterpolatedDistributionForVeloBC(q, f_BN, f_TS, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirTS[k];
+      q = (subgridD.q[TS])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o54* (drho/*+three*(    -vx2+vx3)*/+c9o2*(    -vx2+vx3)*(    -vx2+vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[BN])[kbn]=(c1o1-q)/(c1o1+q)*(f_TS-f_BN+(f_TS+f_BN-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TS+f_BN)-c6o1*c1o54*( -VeloY+VeloZ))/(c1o1+q) - c1o54 * drho;
-         //feq=c1over54* (drho+three*(    -vx2+vx3)+c9over2*(    -vx2+vx3)*(    -vx2+vx3)-cu_sq); 
-         //(D.f[BN])[kbn]=(one-q)/(one+q)*(f_TS-feq*om1)/(one-om1)+(q*(f_TS+f_BN)-six*c1over54*( -VeloY+VeloZ))/(one+q);
-         //(D.f[BN])[kbn]=zero;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+
+         velocityLB = -vx2 + vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o54);
+         velocityBC = -VeloY + VeloZ;
+         (dist.f[BN])[kbn] = getInterpolatedDistributionForVeloBC(q, f_TS, f_BN, feq, omega, velocityBC, c1o54);
       }
 
-      q = q_dirTNE[k];
+      q = (subgridD.q[TNE])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o216*(drho/*+three*( vx1+vx2+vx3)*/+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[BSW])[kbsw]=(c1o1-q)/(c1o1+q)*(f_TNE-f_BSW+(f_TNE+f_BSW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TNE+f_BSW)-c6o1*c1o216*( VeloX+VeloY+VeloZ))/(c1o1+q) - c1o216 * drho;
-         //feq=c1over216*(drho+three*( vx1+vx2+vx3)+c9over2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq); 
-         //(D.f[BSW])[kbsw]=(one-q)/(one+q)*(f_TNE-feq*om1)/(one-om1)+(q*(f_TNE+f_BSW)-six*c1over216*( VeloX+VeloY+VeloZ))/(one+q);
-         //(D.f[BSW])[kbsw]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+         velocityLB = vx1 + vx2 + vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
+         velocityBC = VeloX + VeloY + VeloZ;
+         (dist.f[BSW])[kbsw] = getInterpolatedDistributionForVeloBC(q, f_TNE, f_BSW, feq, omega, velocityBC, c1o216);
       }
 
-      q = q_dirBSW[k];
+      q = (subgridD.q[BSW])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o216*(drho/*+three*(-vx1-vx2-vx3)*/+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[TNE])[ktne]=(c1o1-q)/(c1o1+q)*(f_BSW-f_TNE+(f_BSW+f_TNE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BSW+f_TNE)-c6o1*c1o216*(-VeloX-VeloY-VeloZ))/(c1o1+q) - c1o216 * drho;
-         //feq=c1over216*(drho+three*(-vx1-vx2-vx3)+c9over2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq); 
-         //(D.f[TNE])[ktne]=(one-q)/(one+q)*(f_BSW-feq*om1)/(one-om1)+(q*(f_BSW+f_TNE)-six*c1over216*(-VeloX-VeloY-VeloZ))/(one+q);
-         //(D.f[TNE])[ktne]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+         velocityLB = -vx1 - vx2 - vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
+         velocityBC = -VeloX - VeloY - VeloZ;
+         (dist.f[TNE])[ktne] = getInterpolatedDistributionForVeloBC(q, f_BSW, f_TNE, feq, omega, velocityBC, c1o216);
       }
 
-      q = q_dirBNE[k];
+
+      q = (subgridD.q[BNE])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o216*(drho/*+three*( vx1+vx2-vx3)*/+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[TSW])[ktsw]=(c1o1-q)/(c1o1+q)*(f_BNE-f_TSW+(f_BNE+f_TSW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BNE+f_TSW)-c6o1*c1o216*( VeloX+VeloY-VeloZ))/(c1o1+q) - c1o216 * drho;
-         //feq=c1over216*(drho+three*( vx1+vx2-vx3)+c9over2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq); 
-         //(D.f[TSW])[ktsw]=(one-q)/(one+q)*(f_BNE-feq*om1)/(one-om1)+(q*(f_BNE+f_TSW)-six*c1over216*( VeloX+VeloY-VeloZ))/(one+q);
-         //(D.f[TSW])[ktsw]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+         velocityLB = vx1 + vx2 - vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
+         velocityBC = VeloX + VeloY - VeloZ;
+         (dist.f[TSW])[ktsw] = getInterpolatedDistributionForVeloBC(q, f_BNE, f_TSW, feq, omega, velocityBC, c1o216);
       }
 
-      q = q_dirTSW[k];
+      q = (subgridD.q[TSW])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o216*(drho/*+three*(-vx1-vx2+vx3)*/+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[BNE])[kbne]=(c1o1-q)/(c1o1+q)*(f_TSW-f_BNE+(f_TSW+f_BNE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TSW+f_BNE)-c6o1*c1o216*(-VeloX-VeloY+VeloZ))/(c1o1+q) - c1o216 * drho;
-         //feq=c1over216*(drho+three*(-vx1-vx2+vx3)+c9over2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq); 
-         //(D.f[BNE])[kbne]=(one-q)/(one+q)*(f_TSW-feq*om1)/(one-om1)+(q*(f_TSW+f_BNE)-six*c1over216*(-VeloX-VeloY+VeloZ))/(one+q);
-         //(D.f[BNE])[kbne]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+         velocityLB = -vx1 - vx2 + vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
+         velocityBC = -VeloX - VeloY + VeloZ;
+         (dist.f[BNE])[kbne] = getInterpolatedDistributionForVeloBC(q, f_TSW, f_BNE, feq, omega, velocityBC, c1o216);
       }
 
-      q = q_dirTSE[k];
+      q = (subgridD.q[TSE])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o216*(drho/*+three*( vx1-vx2+vx3)*/+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[BNW])[kbnw]=(c1o1-q)/(c1o1+q)*(f_TSE-f_BNW+(f_TSE+f_BNW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TSE+f_BNW)-c6o1*c1o216*( VeloX-VeloY+VeloZ))/(c1o1+q) - c1o216 * drho;
-         //feq=c1over216*(drho+three*( vx1-vx2+vx3)+c9over2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq); 
-         //(D.f[BNW])[kbnw]=(one-q)/(one+q)*(f_TSE-feq*om1)/(one-om1)+(q*(f_TSE+f_BNW)-six*c1over216*( VeloX-VeloY+VeloZ))/(one+q);
-         //(D.f[BNW])[kbnw]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+         velocityLB = vx1 - vx2 + vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
+         velocityBC = VeloX - VeloY + VeloZ;
+         (dist.f[BNW])[kbnw] = getInterpolatedDistributionForVeloBC(q, f_TSE, f_BNW, feq, omega, velocityBC, c1o216);
       }
 
-      q = q_dirBNW[k];
+      q = (subgridD.q[BNW])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o216*(drho/*+three*(-vx1+vx2-vx3)*/+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[TSE])[ktse]=(c1o1-q)/(c1o1+q)*(f_BNW-f_TSE+(f_BNW+f_TSE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BNW+f_TSE)-c6o1*c1o216*(-VeloX+VeloY-VeloZ))/(c1o1+q) - c1o216 * drho;
-         //feq=c1over216*(drho+three*(-vx1+vx2-vx3)+c9over2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq); 
-         //(D.f[TSE])[ktse]=(one-q)/(one+q)*(f_BNW-feq*om1)/(one-om1)+(q*(f_BNW+f_TSE)-six*c1over216*(-VeloX+VeloY-VeloZ))/(one+q);
-         //(D.f[TSE])[ktse]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+         velocityLB = -vx1 + vx2 - vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
+         velocityBC = -VeloX + VeloY - VeloZ;
+         (dist.f[TSE])[ktse] = getInterpolatedDistributionForVeloBC(q, f_BNW, f_TSE, feq, omega, velocityBC, c1o216);
       }
 
-      q = q_dirBSE[k];
+      q = (subgridD.q[BSE])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o216*(drho/*+three*( vx1-vx2-vx3)*/+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[TNW])[ktnw]=(c1o1-q)/(c1o1+q)*(f_BSE-f_TNW+(f_BSE+f_TNW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BSE+f_TNW)-c6o1*c1o216*( VeloX-VeloY-VeloZ))/(c1o1+q) - c1o216 * drho;
-         //feq=c1over216*(drho+three*( vx1-vx2-vx3)+c9over2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq); 
-         //(D.f[TNW])[ktnw]=(one-q)/(one+q)*(f_BSE-feq*om1)/(one-om1)+(q*(f_BSE+f_TNW)-six*c1over216*( VeloX-VeloY-VeloZ))/(one+q);
-         //(D.f[TNW])[ktnw]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+         velocityLB = vx1 - vx2 - vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
+         velocityBC = VeloX - VeloY - VeloZ;
+         (dist.f[TNW])[ktnw] = getInterpolatedDistributionForVeloBC(q, f_BSE, f_TNW, feq, omega, velocityBC, c1o216);
       }
 
-      q = q_dirTNW[k];
+      q = (subgridD.q[TNW])[k];
       if (q>=c0o1 && q<=c1o1)
       {
-		 VeloX = fac*vx1;
-	     VeloY = fac*vx2;
-	     VeloZ = fac*vx3;
-		 if (x == true) VeloX = c0o1;
-		 if (y == true) VeloY = c0o1;
-		 if (z == true) VeloZ = c0o1;
-         feq=c1o216*(drho/*+three*(-vx1+vx2+vx3)*/+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3) * (c1o1 + drho)-cu_sq); 
-         (D.f[BSE])[kbse]=(c1o1-q)/(c1o1+q)*(f_TNW-f_BSE+(f_TNW+f_BSE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TNW+f_BSE)-c6o1*c1o216*(-VeloX+VeloY+VeloZ))/(c1o1+q) - c1o216 * drho;
-         //feq=c1over216*(drho+three*(-vx1+vx2+vx3)+c9over2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq); 
-         //(D.f[BSE])[kbse]=(one-q)/(one+q)*(f_TNW-feq*om1)/(one-om1)+(q*(f_TNW+f_BSE)-six*c1over216*(-VeloX+VeloY+VeloZ))/(one+q);
-         //(D.f[BSE])[kbse]=zero;
+         VeloX = slipLength*vx1;
+         VeloY = slipLength*vx2;
+         VeloZ = slipLength*vx3;
+         if (x == true) VeloX = c0o1;
+         if (y == true) VeloY = c0o1;
+         if (z == true) VeloZ = c0o1;
+         velocityLB = -vx1 + vx2 + vx3;
+         feq = getEquilibriumForBC(drho, velocityLB, cu_sq, c1o216);
+         velocityBC = -VeloX + VeloY + VeloZ;
+         (dist.f[BSE])[kbse] = getInterpolatedDistributionForVeloBC(q, f_TNW, f_BSE, feq, omega, velocityBC, c1o216);
       }
    }
 }
diff --git a/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.cpp b/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.cpp
index 9c34b930f65335c209e154a472fd3aa31276edf7..334935163712a6e1c304821435b21f2177603b49 100644
--- a/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.cpp
+++ b/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.cpp
@@ -371,40 +371,3 @@ void BCKernelManager::runNoSlipBCKernelPost(const int level) const{
         noSlipBoundaryConditionPost(para->getParD(level).get(), &(para->getParD(level)->noSlipBC));
     }
 }
-
-// void LBKernelManager::calculateMacroscopicValues(const int level) const
-// {
-//     if (para->getIsADcalculationOn()) {
-//           CalcMacADCompSP27(
-//                para->getParD()->velocityX,
-//                para->getParD()->velocityY,
-//                para->getParD()->velocityZ,
-//                para->getParD()->rho,
-//                para->getParD()->pressure,
-//                para->getParD()->typeOfGridNode,
-//                para->getParD()->neighborX,
-//                para->getParD()->neighborY,
-//                para->getParD()->neighborZ,
-//                para->getParD()->numberOfNodes,
-//                para->getParD()->numberofthreads,
-//                para->getParD()->distributions.f[0],
-//                para->getParD()->distributionsAD.f[0],
-//             para->getParD()->forcing,
-//                para->getParD()->isEvenTimestep);
-//     } else {
-//           CalcMacCompSP27(
-//                para->getParD()->velocityX,
-//                para->getParD()->velocityY,
-//                para->getParD()->velocityZ,
-//                para->getParD()->rho,
-//                para->getParD()->pressure,
-//                para->getParD()->typeOfGridNode,
-//                para->getParD()->neighborX,
-//                para->getParD()->neighborY,
-//                para->getParD()->neighborZ,
-//                para->getParD()->numberOfNodes,
-//                para->getParD()->numberofthreads,
-//                para->getParD()->distributions.f[0],
-//                para->getParD()->isEvenTimestep);
-//      }
-// }
diff --git a/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.h b/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.h
index d4d87af7612fccf94fd2ae5e4c92cfcac01c6013..bc4cbb056abe86b817d6b784135c930fa89660c4 100644
--- a/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.h
+++ b/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.h
@@ -26,12 +26,12 @@
 //  You should have received a copy of the GNU General Public License along
 //  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
 //
-//! \file LBKernelManager.h
+//! \file BCKernelManager.h
 //! \ingroup KernelManager
 //! \author Martin Schoenherr
 //=======================================================================================
-#ifndef LBKernelManager_H
-#define LBKernelManager_H
+#ifndef BCKernelManager_H
+#define BCKernelManager_H
 
 #include <functional>
 #include <memory>
@@ -48,8 +48,8 @@ struct LBMSimulationParameter;
 using boundaryCondition = std::function<void(LBMSimulationParameter *, QforBoundaryConditions *)>;
 using boundaryConditionPara = std::function<void(Parameter *, QforBoundaryConditions *, const int level)>;
 
-//! \class LBKernelManager
-//! \brief manage the cuda kernel calls
+//! \class BCKernelManager
+//! \brief manage the cuda kernel calls to boundary conditions
 class VIRTUALFLUIDS_GPU_EXPORT BCKernelManager
 {
 public:
@@ -57,8 +57,6 @@ public:
     //! \param parameter shared pointer to instance of class Parameter
     BCKernelManager(SPtr<Parameter> parameter, BoundaryConditionFactory *bcFactory);
 
-    void setBoundaryConditionKernels();
-
     //! \brief calls the device function of the velocity boundary condition (post-collision)
     void runVelocityBCKernelPost(const int level) const;
 
@@ -71,10 +69,10 @@ public:
     //! \brief calls the device function of the geometry boundary condition (pre-collision)
     void runGeoBCKernelPre(const int level, unsigned int t, CudaMemoryManager *cudaMemoryManager) const;
 
-    //! \brief calls the device function of the slip boundary condition
+    //! \brief calls the device function of the slip boundary condition (post-collision)
     void runSlipBCKernelPost(const int level) const;
 
-    //! \brief calls the device function of the no-slip boundary condition
+    //! \brief calls the device function of the no-slip boundary condition (post-collision)
     void runNoSlipBCKernelPost(const int level) const;
 
     //! \brief calls the device function of the pressure boundary condition (pre-collision)
@@ -83,15 +81,12 @@ public:
     //! \brief calls the device function of the pressure boundary condition (post-collision)
     void runPressureBCKernelPost(const int level) const;
 
-    //! \brief calls the device function of the outflow boundary condition
+    //! \brief calls the device function of the outflow boundary condition (pre-collision)
     void runOutflowBCKernelPre(const int level) const;
 
-    //! \brief calls the device function of the stress wall model
+    //! \brief calls the device function of the stress wall model (post-collision)
     void runStressWallModelKernelPost(const int level) const;
 
-    //! \brief calls the device function that calculates the macroscopic values
-    void calculateMacroscopicValues(const int level) const;
-
 private:
     SPtr<Parameter> para;
 
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index 4cdb9f1db6514d61ff90327a928822df924e7270..a22cd110e9b29019e5c52aab256cb5bb2a0103fe 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -1,6 +1,5 @@
 #include "Simulation.h"
 
-#include <stdio.h>
 #include <vector>
 
 #include <helper_timer.h>
@@ -44,6 +43,7 @@
 #include "Calculation/PorousMedia.h"
 //////////////////////////////////////////////////////////////////////////
 #include "Output/Timer.h"
+#include "Output/FileWriter.h"
 //////////////////////////////////////////////////////////////////////////
 #include "Restart/RestartObject.h"
 //////////////////////////////////////////////////////////////////////////
@@ -59,7 +59,6 @@
 
 #include <logger/Logger.h>
 
-#include "Output/FileWriter.h"
 
 
 std::string getFileName(const std::string& fname, int step, int myID)
@@ -173,7 +172,7 @@ Simulation::Simulation(std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemo
     //////////////////////////////////////////////////////////////////////////
     // Particles preprocessing
     //////////////////////////////////////////////////////////////////////////
-    if (para->getCalcParticle()) {
+    if (para->getCalcParticles()) {
         rearrangeGeometry(para.get(), cudaMemoryManager.get());
         //////////////////////////////////////////////////////////////////////////
         allocParticles(para.get(), cudaMemoryManager.get());
@@ -384,7 +383,7 @@ Simulation::Simulation(std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemo
     //////////////////////////////////////////////////////////////////////////
     output << "Print files Init...";
     dataWriter->writeInit(para, cudaMemoryManager);
-    if (para->getCalcParticle())
+    if (para->getCalcParticles())
         copyAndPrintParticles(para.get(), cudaMemoryManager.get(), 0, true);
     output << "done.\n";
 
@@ -442,7 +441,7 @@ void Simulation::allocNeighborsOffsetsScalesAndBoundaries(GridProvider &gridProv
 
 void Simulation::run()
 {
-   unsigned int t, t_prev;
+   unsigned int timestep, t_prev;
    uint t_turbulenceIntensity = 0;
    unsigned int t_MP = 0;
 
@@ -473,14 +472,14 @@ void Simulation::run()
 	////////////////////////////////////////////////////////////////////////////////
 	// Time loop
 	////////////////////////////////////////////////////////////////////////////////
-	for(t=para->getTStart();t<=para->getTEnd();t++)
+	for(timestep=para->getTStart();timestep<=para->getTEnd();timestep++)
 	{
-        this->updateGrid27->updateGrid(0, t);
+        this->updateGrid27->updateGrid(0, timestep);
 
 	    ////////////////////////////////////////////////////////////////////////////////
 	    //Particles
 	    ////////////////////////////////////////////////////////////////////////////////
-	    if (para->getCalcParticle()) propagateParticles(para.get(), t);
+	    if (para->getCalcParticles()) propagateParticles(para.get(), timestep);
 	    ////////////////////////////////////////////////////////////////////////////////
 
 
@@ -495,8 +494,8 @@ void Simulation::run()
             exchangeMultiGPU(para.get(), communicator, cudaMemoryManager.get(), 0, -1);
         }
 
-	    if( this->kineticEnergyAnalyzer ) this->kineticEnergyAnalyzer->run(t);
-	    if( this->enstrophyAnalyzer     ) this->enstrophyAnalyzer->run(t);
+	    if( this->kineticEnergyAnalyzer ) this->kineticEnergyAnalyzer->run(timestep);
+	    if( this->enstrophyAnalyzer     ) this->enstrophyAnalyzer->run(timestep);
 	    ////////////////////////////////////////////////////////////////////////////////
 
 
@@ -505,7 +504,7 @@ void Simulation::run()
         ////////////////////////////////////////////////////////////////////////////////
         //Calc Median
         ////////////////////////////////////////////////////////////////////////////////
-        if (para->getCalcMedian() && ((int)t >= para->getTimeCalcMedStart()) && ((int)t <= para->getTimeCalcMedEnd()))
+        if (para->getCalcMedian() && ((int)timestep >= para->getTimeCalcMedStart()) && ((int)timestep <= para->getTimeCalcMedEnd()))
         {
           for (int lev=para->getCoarse(); lev <= para->getFine(); lev++)
           {
@@ -573,23 +572,23 @@ void Simulation::run()
         ////////////////////////////////////////////////////////////////////////////////
         // CheckPoint
         ////////////////////////////////////////////////////////////////////////////////
-        if(para->getDoCheckPoint() && para->getTimeDoCheckPoint()>0 && t%para->getTimeDoCheckPoint()==0 && t>0 && !para->overWritingRestart(t))
+        if(para->getDoCheckPoint() && para->getTimeDoCheckPoint()>0 && timestep%para->getTimeDoCheckPoint()==0 && timestep>0 && !para->overWritingRestart(timestep))
         {
 			averageTimer->stopTimer();
             //////////////////////////////////////////////////////////////////////////
 
             if( para->getDoCheckPoint() )
             {
-                output << "Copy data for CheckPoint t=" << t << "...\n";
+                output << "Copy data for CheckPoint t=" << timestep << "...\n";
 
                 for (int lev=para->getCoarse(); lev <= para->getFine(); lev++)
                 {
                     cudaMemoryManager->cudaCopyFsForCheckPoint(lev);
                 }
 
-                output << "Write data for CheckPoint t=" << t << "...";
+                output << "Write data for CheckPoint t=" << timestep << "...";
 
-				const auto name = getFileName(para->getFName(), t, para->getMyID());
+				const auto name = getFileName(para->getFName(), timestep, para->getMyID());
 				restart_object->serialize(name, para);
 
                 output << "\n done\n";
@@ -609,7 +608,7 @@ void Simulation::run()
         //set MP-Time
         if (para->getUseMeasurePoints())
         {
-            if ((t%para->getTimestepForMP()) == 0)
+            if ((timestep%para->getTimestepForMP()) == 0)
             {
                 unsigned int valuesPerClockCycle = (unsigned int)(para->getclockCycleForMP() / para->getTimestepForMP());
                 for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
@@ -626,16 +625,16 @@ void Simulation::run()
             }
 
             //Copy Measure Values
-            if ((t % (unsigned int)para->getclockCycleForMP()) == 0)
+            if ((timestep % (unsigned int)para->getclockCycleForMP()) == 0)
             {
                 for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
                 {
                     cudaMemoryManager->cudaCopyMeasurePointsToHost(lev);
                     para->copyMeasurePointsArrayToVector(lev);
-                    output << "\n Write MeasurePoints at level = " << lev << " and timestep = " << t << "\n";
+                    output << "\n Write MeasurePoints at level = " << lev << " and timestep = " << timestep << "\n";
                     for (int j = 0; j < (int)para->getParH(lev)->MP.size(); j++)
                     {
-                        MeasurePointWriter::writeMeasurePoints(para.get(), lev, j, t);
+                        MeasurePointWriter::writeMeasurePoints(para.get(), lev, j, timestep);
                     }
                     //MeasurePointWriter::calcAndWriteMeanAndFluctuations(para.get(), lev, t, para->getTStartOut());
                 }
@@ -702,7 +701,7 @@ void Simulation::run()
       // File IO
       ////////////////////////////////////////////////////////////////////////////////
       //communicator->startTimer();
-      if(para->getTOut()>0 && t%para->getTOut()==0 && t>para->getTStartOut())
+      if(para->getTOut()>0 && timestep%para->getTOut()==0 && timestep>para->getTStartOut())
       {
 		  //////////////////////////////////////////////////////////////////////////////////
 		  //if (para->getParD(0)->evenOrOdd==true)  para->getParD(0)->evenOrOdd=false;
@@ -711,12 +710,12 @@ void Simulation::run()
 
 		//////////////////////////////////////////////////////////////////////////
 		averageTimer->stopTimer();
-		averageTimer->outputPerformance(t, para.get(), communicator);
+		averageTimer->outputPerformance(timestep, para.get(), communicator);
 		//////////////////////////////////////////////////////////////////////////
 
          if( para->getPrintFiles() )
          {
-            output << "Write files t=" << t << "... ";
+            output << "Write files t=" << timestep << "... ";
             for (int lev=para->getCoarse(); lev <= para->getFine(); lev++)
             {
 		        //////////////////////////////////////////////////////////////////////////
@@ -834,11 +833,11 @@ void Simulation::run()
 			   //////////////////////////////////////////////////////////////////////////
                //TODO: implement flag to write ASCII data
 			   if (para->getWriteVeloASCIIfiles())
-				   VeloASCIIWriter::writeVelocitiesAsTXT(para.get(), lev, t);
+				   VeloASCIIWriter::writeVelocitiesAsTXT(para.get(), lev, timestep);
 			   //////////////////////////////////////////////////////////////////////////
                if( this->kineticEnergyAnalyzer || this->enstrophyAnalyzer )
                {
-                   std::string fname = para->getFName() + "_ID_" + StringUtil::toString<int>(para->getMyID()) + "_t_" + StringUtil::toString<int>(t);
+                   std::string fname = para->getFName() + "_ID_" + StringUtil::toString<int>(para->getMyID()) + "_t_" + StringUtil::toString<int>(timestep);
 
                    if (this->kineticEnergyAnalyzer) this->kineticEnergyAnalyzer->writeToFile(fname);
                    if (this->enstrophyAnalyzer)     this->enstrophyAnalyzer->writeToFile(fname);
@@ -964,34 +963,34 @@ void Simulation::run()
 			////////////////////////////////////////////////////////////////////////
 			//calculate median on host
 			////////////////////////////////////////////////////////////////////////
-			if (para->getCalcMedian() && ((int)t > para->getTimeCalcMedStart()) && ((int)t <= para->getTimeCalcMedEnd()) && ((t%(unsigned int)para->getclockCycleForMP())==0))
+			if (para->getCalcMedian() && ((int)timestep > para->getTimeCalcMedStart()) && ((int)timestep <= para->getTimeCalcMedEnd()) && ((timestep%(unsigned int)para->getclockCycleForMP())==0))
 			{
-				unsigned int tdiff = t - t_prev;
+				unsigned int tdiff = timestep - t_prev;
 				calcMedian(para.get(), tdiff);
 
 				/////////////////////////////////
 				//added for incremental averaging
-				t_prev = t;
+				t_prev = timestep;
 				resetMedian(para.get());
 				/////////////////////////////////
 			}
             if (para->getCalcTurbulenceIntensity())
 			{
-                uint t_diff = t - t_turbulenceIntensity;
+                uint t_diff = timestep - t_turbulenceIntensity;
                 calcTurbulenceIntensity(para.get(), cudaMemoryManager.get(), t_diff);
                 //writeAllTiDatafToFile(para.get(), t);
             }
 			////////////////////////////////////////////////////////////////////////
-			dataWriter->writeTimestep(para, t);
+			dataWriter->writeTimestep(para, timestep);
 			////////////////////////////////////////////////////////////////////////
             if (para->getCalcTurbulenceIntensity()) {
-                t_turbulenceIntensity = t;
+                t_turbulenceIntensity = timestep;
                 resetVelocityFluctuationsAndMeans(para.get(), cudaMemoryManager.get());
             }
 			////////////////////////////////////////////////////////////////////////
-            if (para->getCalcDragLift()) printDragLift(para.get(), cudaMemoryManager.get(), t);
+            if (para->getCalcDragLift()) printDragLift(para.get(), cudaMemoryManager.get(), timestep);
 			////////////////////////////////////////////////////////////////////////
-			if (para->getCalcParticle()) copyAndPrintParticles(para.get(), cudaMemoryManager.get(), t, false);
+			if (para->getCalcParticles()) copyAndPrintParticles(para.get(), cudaMemoryManager.get(), timestep, false);
 			////////////////////////////////////////////////////////////////////////
 			output << "done.\n";
 			////////////////////////////////////////////////////////////////////////
diff --git a/src/gpu/VirtualFluids_GPU/Output/Timer.cpp b/src/gpu/VirtualFluids_GPU/Output/Timer.cpp
index 3f95d95a611da3c75e57dcb3431024e1ca71bd27..0a8ef510d3b1ec2137d56538f55423f9940c2037 100644
--- a/src/gpu/VirtualFluids_GPU/Output/Timer.cpp
+++ b/src/gpu/VirtualFluids_GPU/Output/Timer.cpp
@@ -1,11 +1,11 @@
-
+#include "Timer.h"
 #include <iostream>
-#include <cuda_runtime.h>
+#include <helper_cuda.h>
+
 #include "UbScheduler.h"
-#include "Timer.h"
+#include "Parameter/Parameter.h"
 #include "VirtualFluids_GPU/Communication/Communicator.h"
 
-
 void Timer::initTimer()
 {
     cudaEventCreate(&this->start_t);
diff --git a/src/gpu/VirtualFluids_GPU/Output/Timer.h b/src/gpu/VirtualFluids_GPU/Output/Timer.h
index 26be785c7f76b7695656c9600bdb586804dca251..d035cbb6cef7ea9f8edabbd2894671a868c37eec 100644
--- a/src/gpu/VirtualFluids_GPU/Output/Timer.h
+++ b/src/gpu/VirtualFluids_GPU/Output/Timer.h
@@ -1,17 +1,15 @@
 #ifndef TIMER_H
 #define TIMER_H
-
-#include "helper_cuda.h"
 #include <cuda_runtime.h>
-#include "Core/DataTypes.h"
 
-#include "UbScheduler.h"
-#include "logger/Logger.h"
+#include "Core/DataTypes.h"
 #include "Parameter/Parameter.h"
+#include "logger/Logger.h"
 
 namespace vf::gpu{
     class Communicator;
 }
+class Parameter;
 
 class Timer
 {
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/EdgeNodeFinder.cpp b/src/gpu/VirtualFluids_GPU/Parameter/EdgeNodeFinder.cpp
index a2662249f25b30dda9d582e188d5ff126b2c2c5b..fc1d9b19752be009c8105e5d2759c33afe8e9400 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/EdgeNodeFinder.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/EdgeNodeFinder.cpp
@@ -7,11 +7,11 @@
 namespace vf::gpu
 {
 //! \brief Find nodes that are both received in the x-direction and sent in the y-direction
-void findEdgeNodesXY(const int level, LBMSimulationParameter& parameterLB);
+void findEdgeNodesXY(LBMSimulationParameter& parameterLB);
 //! \brief Find nodes that are both received in the x-direction and sent in the z-direction
-void findEdgeNodesXZ(const int level, LBMSimulationParameter& parameterLB);
+void findEdgeNodesXZ(LBMSimulationParameter& parameterLB);
 //! \brief Find nodes that are both received in the y-direction and sent in the z-direction
-void findEdgeNodesYZ(const int level, LBMSimulationParameter& parameterLB);
+void findEdgeNodesYZ(LBMSimulationParameter& parameterLB);
 void findEdgeNodes(const std::vector<ProcessNeighbor27> &recvProcessNeighbor,
                    const std::vector<ProcessNeighbor27> &sendProcessNeighbor,
                    std::vector<LBMSimulationParameter::EdgeNodePositions> &edgeNodes);
@@ -21,25 +21,25 @@ std::optional<std::pair<int, int>> findIndexInSendNodes(const int nodeIndex,
 void findEdgeNodesCommMultiGPU(Parameter& parameter)
 {
     for (int level = 0; level <= parameter.getFine(); level++) {
-        findEdgeNodesXY(level, *parameter.getParH(level));
-        findEdgeNodesXZ(level, *parameter.getParH(level));
-        findEdgeNodesYZ(level, *parameter.getParH(level));
+        findEdgeNodesXY(*parameter.getParH(level));
+        findEdgeNodesXZ(*parameter.getParH(level));
+        findEdgeNodesYZ(*parameter.getParH(level));
     }
 }
 
-void findEdgeNodesXY(const int level, LBMSimulationParameter& parameterLB)
+void findEdgeNodesXY(LBMSimulationParameter& parameterLB)
 {
     findEdgeNodes(parameterLB.recvProcessNeighborX, parameterLB.sendProcessNeighborY,
                   parameterLB.edgeNodesXtoY);
 }
 
-void findEdgeNodesXZ(const int level, LBMSimulationParameter& parameterLB)
+void findEdgeNodesXZ(LBMSimulationParameter& parameterLB)
 {
     findEdgeNodes(parameterLB.recvProcessNeighborX, parameterLB.sendProcessNeighborZ,
                   parameterLB.edgeNodesXtoZ);
 }
 
-void findEdgeNodesYZ(const int level, LBMSimulationParameter& parameterLB)
+void findEdgeNodesYZ(LBMSimulationParameter& parameterLB)
 {
     findEdgeNodes(parameterLB.recvProcessNeighborY, parameterLB.sendProcessNeighborZ,
                   parameterLB.edgeNodesYtoZ);
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
index 42e53bb681b3977f42ea8ebfebd0fbdebae37444..5ae009f589c76b5307a653b9b1b9e385e2d848ae 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -92,7 +92,7 @@ void Parameter::readConfigData(const vf::basics::ConfigurationFile &configData)
         this->setPrintFiles(configData.getValue<bool>("WriteGrid"));
     //////////////////////////////////////////////////////////////////////////
     if (configData.contains("GeometryValues"))
-        this->setGeometryValues(configData.getValue<bool>("GeometryValues"));
+        this->setUseGeometryValues(configData.getValue<bool>("GeometryValues"));
     //////////////////////////////////////////////////////////////////////////
     if (configData.contains("calc2ndOrderMoments"))
         this->setCalc2ndOrderMoments(configData.getValue<bool>("calc2ndOrderMoments"));
@@ -146,10 +146,10 @@ void Parameter::readConfigData(const vf::basics::ConfigurationFile &configData)
         this->setD3Qxx(configData.getValue<int>("D3Qxx"));
     //////////////////////////////////////////////////////////////////////////
     if (configData.contains("TimeEnd"))
-        this->setTEnd(configData.getValue<int>("TimeEnd"));
+        this->setTimestepEnd(configData.getValue<int>("TimeEnd"));
     //////////////////////////////////////////////////////////////////////////
     if (configData.contains("TimeOut"))
-        this->setTOut(configData.getValue<int>("TimeOut"));
+        this->setTimestepOut(configData.getValue<int>("TimeOut"));
     //////////////////////////////////////////////////////////////////////////
     if (configData.contains("TimeStartOut"))
         this->setTStartOut(configData.getValue<int>("TimeStartOut"));
@@ -191,10 +191,10 @@ void Parameter::readConfigData(const vf::basics::ConfigurationFile &configData)
 
     //////////////////////////////////////////////////////////////////////////
     if (configData.contains("Viscosity_LB"))
-        this->setViscosity(configData.getValue<real>("Viscosity_LB"));
+        this->setViscosityLB(configData.getValue<real>("Viscosity_LB"));
     //////////////////////////////////////////////////////////////////////////
     if (configData.contains("Velocity_LB"))
-        this->setVelocity(configData.getValue<real>("Velocity_LB"));
+        this->setVelocityLB(configData.getValue<real>("Velocity_LB"));
     //////////////////////////////////////////////////////////////////////////
     if (configData.contains("Viscosity_Ratio_World_to_LB"))
         this->setViscosityRatio(configData.getValue<real>("Viscosity_Ratio_World_to_LB"));
@@ -699,11 +699,11 @@ void Parameter::setEndXHotWall(real endXHotWall)
 {
     this->endXHotWall = endXHotWall;
 }
-void Parameter::setTEnd(unsigned int tend)
+void Parameter::setTimestepEnd(unsigned int tend)
 {
     ic.tend = tend;
 }
-void Parameter::setTOut(unsigned int tout)
+void Parameter::setTimestepOut(unsigned int tout)
 {
     ic.tout = tout;
 }
@@ -788,11 +788,11 @@ void Parameter::setTemperatureBC(real TempBC)
 {
     ic.TempBC = TempBC;
 }
-void Parameter::setViscosity(real Viscosity)
+void Parameter::setViscosityLB(real Viscosity)
 {
     ic.vis = Viscosity;
 }
-void Parameter::setVelocity(real Velocity)
+void Parameter::setVelocityLB(real Velocity)
 {
     ic.u0 = Velocity;
 }
@@ -1331,9 +1331,9 @@ void Parameter::setObj(std::string str, bool isObj)
         this->setIsOutflowNormal(isObj);
     }
 }
-void Parameter::setGeometryValues(bool GeometryValues)
+void Parameter::setUseGeometryValues(bool useGeometryValues)
 {
-    ic.GeometryValues = GeometryValues;
+    ic.GeometryValues = useGeometryValues;
 }
 void Parameter::setCalc2ndOrderMoments(bool is2ndOrderMoments)
 {
@@ -1781,7 +1781,7 @@ bool Parameter::getCalcCp()
 {
     return this->calcCp;
 }
-bool Parameter::getCalcParticle()
+bool Parameter::getCalcParticles()
 {
     return this->calcParticles;
 }
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
index c6adcac6a9ae83308d843d3bdc819f68f0e4ca03..c1386df80e80082e9a31faeaa593d0c0aca7d06a 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
@@ -374,7 +374,7 @@ struct LBMSimulationParameter {
 class VIRTUALFLUIDS_GPU_EXPORT Parameter
 {
 public:
-    Parameter(const vf::basics::ConfigurationFile &configData, const int numberOfProcesses = 1, const int myId = 0);
+    explicit Parameter(const vf::basics::ConfigurationFile &configData, const int numberOfProcesses = 1, const int myId = 0);
     Parameter(const int numberOfProcesses = 1, const int myId = 0);
     ~Parameter();
     void initLBMSimulationParameter();
@@ -408,8 +408,8 @@ public:
     void setCalcParticles(bool calcParticles);
     void setStartXHotWall(real startXHotWall);
     void setEndXHotWall(real endXHotWall);
-    void setTEnd(unsigned int tend);
-    void setTOut(unsigned int tout);
+    void setTimestepEnd(unsigned int tend);
+    void setTimestepOut(unsigned int tout);
     void setTStartOut(unsigned int tStartOut);
     void setTimestepOfCoarseLevel(unsigned int timestep);
     void setCalcTurbulenceIntensity(bool calcVelocityAndFluctuations);
@@ -497,8 +497,8 @@ public:
     void setReadGeo(bool readGeo);
     void setTemperatureInit(real Temp);
     void setTemperatureBC(real TempBC);
-    void setViscosity(real Viscosity);
-    void setVelocity(real Velocity);
+    void setViscosityLB(real Viscosity);
+    void setVelocityLB(real Velocity);
     void setViscosityRatio(real ViscosityRatio);
     void setVelocityRatio(real VelocityRatio);
     void setDensityRatio(real DensityRatio);
@@ -552,7 +552,7 @@ public:
     void setDoCheckPoint(bool doCheckPoint);
     void setDoRestart(bool doRestart);
     void setObj(std::string str, bool isObj);
-    void setGeometryValues(bool GeometryValues);
+    void setUseGeometryValues(bool GeometryValues);
     void setCalc2ndOrderMoments(bool is2ndOrderMoments);
     void setCalc3rdOrderMoments(bool is3rdOrderMoments);
     void setCalcHighOrderMoments(bool isHighOrderMoments);
@@ -627,7 +627,7 @@ public:
     bool getCalcMedian();
     bool getCalcDragLift();
     bool getCalcCp();
-    bool getCalcParticle();
+    bool getCalcParticles();
     bool getWriteVeloASCIIfiles();
     bool getCalcPlaneConc();
     //! \returns index of finest level
@@ -883,7 +883,7 @@ private:
     unsigned int timestep;
 
     // Kernel
-    std::string mainKernel{ "CumulantK17Comp" };
+    std::string mainKernel{ "CumulantK17CompChim" };
     bool multiKernelOn{ false };
     std::vector<int> multiKernelLevel;
     std::vector<std::string> multiKernel;
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/ParameterTest.cpp b/src/gpu/VirtualFluids_GPU/Parameter/ParameterTest.cpp
index 9e05ed1332b34420656e6c0c81f07501da7c7aac..b2d1306e2ab262224cf2b000ccc792618fd98c5a 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/ParameterTest.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/ParameterTest.cpp
@@ -113,7 +113,7 @@ TEST(ParameterTest, check_all_Parameter_CanBePassedToConstructor)
         EXPECT_THAT((real)limiters_actual[i], RealEq(limiters[i]));
     }
 
-    EXPECT_THAT(para.getCalcParticle(), testing::Eq(true));
+    EXPECT_THAT(para.getCalcParticles(), testing::Eq(true));
     EXPECT_THAT(para.getParticleBasicLevel(), testing::Eq(1));
     EXPECT_THAT(para.getParticleInitLevel(), testing::Eq(2));
     EXPECT_THAT(para.getNumberOfParticles(), testing::Eq(1111));