diff --git a/CMake/cmake_config_files/MOLLOK.config.cmake b/CMake/cmake_config_files/MOLLOK.config.cmake
index 2407cda6287cd0e171453d1559ab6585b0a09a33..03f83455175719ff5b0e786994213bbaa0fbd29e 100644
--- a/CMake/cmake_config_files/MOLLOK.config.cmake
+++ b/CMake/cmake_config_files/MOLLOK.config.cmake
@@ -5,7 +5,7 @@
 #################################################################################
 
 #SET TO CORRECT PATH:
-SET(CMAKE_CUDA_ARCHITECTURES 52)
+SET(CMAKE_CUDA_ARCHITECTURES 86)
 
 SET(PATH_NUMERICAL_TESTS "D:/out/numericalTests/")
 LIST(APPEND VF_COMPILER_DEFINITION "PATH_NUMERICAL_TESTS=${PATH_NUMERICAL_TESTS}")
diff --git a/CMake/cmake_config_files/TESLA03.config.cmake b/CMake/cmake_config_files/TESLA03.config.cmake
index bf08ef36405df6090393f368e50e4e48795635af..91672511794eeb452f97a1aca952416c7e1179d5 100644
--- a/CMake/cmake_config_files/TESLA03.config.cmake
+++ b/CMake/cmake_config_files/TESLA03.config.cmake
@@ -19,3 +19,4 @@ SET(VTK_DIR "F:/Libraries/vtk/VTK-8.2.0/build" CACHE PATH "VTK directory overrid
 
 SET(PATH_NUMERICAL_TESTS "E:/temp/numericalTests/")
 LIST(APPEND VF_COMPILER_DEFINITION "PATH_NUMERICAL_TESTS=${PATH_NUMERICAL_TESTS}")
+SET(CMAKE_CUDA_ARCHITECTURES 52)
\ No newline at end of file
diff --git a/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp b/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp
index 2263bf9ffec4c4e215d9abd622104040f31fa707..86d2227d200622ab3a2d87505ba5c7afcede36fe 100644
--- a/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp
+++ b/apps/gpu/LBM/WTG_RUB/WTG_RUB.cpp
@@ -94,6 +94,8 @@ const uint timeStepEnd = 100000;
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+void convertMPFile(SPtr <MultipleGridBuilder> gridBuilder, real& phi, std::vector<real>& origin, bool& measureVeloProfilesOnly, uint& maxLevel);
+
 void addFineGrids(SPtr<MultipleGridBuilder> gridBuilder, uint &maxLevel, real &rotationOfCity);
 
 void readVelocityProfile();
@@ -143,7 +145,7 @@ void multipleLevel(const std::string& configPath)
 
     //TriangularMesh *RubSTL      = TriangularMesh::make(inputPath + "stl/Var02_0deg_FD_b.stl");
     TriangularMesh *RubSTL      = TriangularMesh::make(inputPath + "stl/" + chooseVariation() + ".stl");
-    // vector<real> originOfCityXY = { 600.0, y_max / 2, z_offset };
+    std::vector<real> originOfCityXY = { 600.0, y_max / 2, z_offset };
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     // OPTIONS
@@ -210,7 +212,7 @@ void multipleLevel(const std::string& configPath)
 	para->setDevices(std::vector<uint>{(uint)0});
 
     para->setOutputPath( path );
-    para->setOutputPrefix( simulationName );
+    para->setOutputPrefix( "Unified_" + simulationName );
 
     para->setFName(para->getOutputPath() + "/" + para->getOutputPrefix());
 
@@ -223,7 +225,7 @@ void multipleLevel(const std::string& configPath)
 
     para->setVelocityRatio(velocity/ velocityLB);
 
-	para->setMainKernel("CumulantK17CompChim");
+	para->setMainKernel("CumulantK17CompChim"); // CumulantK17Unified, CumulantK17CompChim
 
 	para->setInitialCondition([&](real coordX, real coordY, real coordZ, real &rho, real &vx, real &vy, real &vz) {
         rho = (real)0.0;
@@ -305,6 +307,18 @@ void multipleLevel(const std::string& configPath)
         });
     }
 
+    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Intializing MeasurePoints
+    para->setUseMeasurePoints(useMP);
+    if (para->getUseMeasurePoints()) {
+        convertMPFile(gridBuilder, rotationOfCity, originOfCityXY, measureVeloProfilesOnly, maxLevel);
+        // Number log-Files for each MeasurePoint: numberOfMPFiles = timeStepEnd/ClockCycle
+        para->setclockCycleForMP(timeStepEnd);
+        // Number of  logged Timesteps for each file
+        para->settimestepForMP(timeStepOut / 100);
+        // para->settimestepForMP(timeStepOut);
+        para->setmeasurePoints(inputPath + "measurePoints.dat");
+    }
 
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -522,6 +536,131 @@ void addFineGrids(SPtr<MultipleGridBuilder> gridBuilder, uint &maxLevel, real &r
 
 }
 
+void convertMPFile(SPtr<MultipleGridBuilder> gridBuilder, real &phi, std::vector<real> &originXY, bool &measureVeloProfilesOnly, uint &maxLevel)
+{
+    // File Reader&Writer for converting MP-Coordinates to Index: MeasurePoint placement requires "measurePoints.dat"
+    // with name, node-ID and level. This function can read a txt-File containing the name, X-Y-Z-Coordinates and level
+    // of measurePoints. After reading the txt-File and converting X-Y-Z to the node-ID, it writes "measurePoints.dat".
+    // Justification for this function: Human Readability and no changes in measurepoint core functions
+
+    // File Opening Procedure
+    std::ifstream inFile;
+    if (measureVeloProfilesOnly)
+        inFile.open(inputPath + "measurePoints_veloProfiles.txt");
+    else
+        inFile.open(inputPath + "measurePoints.txt");
+
+    // Check for error
+    if (inFile.fail()) {
+        std::cerr << "Error opening File" << std::endl;
+        exit(1);
+    }
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Reading Procedure
+    std::cout << "phi in degrees:" << phi << std::endl;
+    phi = phi * M_PI / 180;
+    std::cout << "phi in radians:" << phi << std::endl;
+
+    std::vector<std::string> MP_name;
+    std::vector<real> MP_X, MP_Y, MP_Z;
+    std::vector<int> MP_level, MP_k;
+
+    std::string name;
+    real X, Y, Z;
+    uint level, numberOfMeasurePoints;
+
+    inFile >> numberOfMeasurePoints;
+    std::cout << "numberOfMeasurePoints: " << numberOfMeasurePoints << " ";
+    std::cout << "Coordinates from File\n";
+    for (uint k = 0; k < numberOfMeasurePoints; k++) {
+        inFile >> name;
+        MP_name.push_back(name);
+        std::cout << "Name: " << MP_name[k] << " ";
+
+        inFile >> X;
+        MP_X.push_back(X);
+        std::cout << "\t\tX: " << MP_X[k] << " ";
+
+        inFile >> Y;
+        MP_Y.push_back(Y);
+        std::cout << "\t\tY: " << MP_Y[k] << " ";
+
+        inFile >> Z;
+        if (((variant > 3 && variant < 7) || (variant > 9 && variant <= 12)) && k == 14)
+            Z += 2.25; // account for angled roof
+        MP_Z.push_back(Z);
+        std::cout << "\t\tZ: " << MP_Z[k] + z_offset << " ";
+
+        inFile >> level;
+        if (level > maxLevel)
+            level = maxLevel;
+        MP_level.push_back(level);
+        std::cout << "\t\tLevel: " << MP_level[k] << std::endl;
+    }
+    inFile.close();
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    real X_temp, Y_temp;
+    // Transformation for phi radians around centre of city for MP[1..15]
+    if (!phi == 0) {
+        std::cout << "Calculating new Coordinates for MP01 to MP15 after Rotation of " << phi * 180 / M_PI
+                  << "degrees (+: counter-clockwise / -: clockwise)\n";
+        for (uint k = 0; k < 15; k++) {
+            X_temp = originXY[0] + (MP_X[k] - originXY[0]) * cos(phi) - (MP_Y[k] - originXY[1]) * sin(phi);
+            Y_temp = originXY[1] + (MP_X[k] - originXY[0]) * sin(phi) + (MP_Y[k] - originXY[1]) * cos(phi);
+            std::cout << "Name:  " << MP_name[k] << " ";
+            std::cout << "\t\tX: " << X_temp << " ";
+            std::cout << "\t\tY: " << Y_temp << " ";
+            std::cout << "\t\tZ: " << MP_Z[k] << " " << std::endl;
+
+            MP_X[k] = X_temp;
+            MP_Y[k] = Y_temp;
+        }
+    }
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Coordinates to Index Procedure
+    // std::cout << "Converting Coordinates to Index..." << std::endl;
+    for (uint k = 0; k < numberOfMeasurePoints; k++) {
+        MP_k.push_back(
+            gridBuilder->getGrid(MP_level[k])
+                ->getSparseIndex(gridBuilder->getGrid(MP_level[k])->transCoordToIndex(MP_X[k], MP_Y[k], MP_Z[k])));
+        if (MP_k[k] == -1) {
+            std::cerr << "Error: Could not convert Coordinate to Sparse Index for MP " << k + 1 << std::endl;
+        }
+        std::cout << MP_name[k] << "\tID = "
+                  << gridBuilder->getGrid(MP_level[k])
+                         ->getSparseIndex(
+                             gridBuilder->getGrid(MP_level[k])->transCoordToIndex(MP_X[k], MP_Y[k], MP_Z[k]))
+                  << std::endl;
+        // std::cout << "ID = " <<
+        // gridBuilder->getGrid(0)->getSparseIndex(gridBuilder->getGrid(0)->transCoordToIndex(-0.500000,
+        // -0.500000, 9.500000)) << std::endl;
+    }
+    // std::cout << "Done Converting Coordinates to Index..." << std::endl;
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Writing Procedure
+    // std::cout << "Writing new file..." << std::endl;
+    std::ofstream outFile(inputPath + "measurePoints.dat");
+
+    outFile << numberOfMeasurePoints << std::endl;
+    for (uint j = 0; j < numberOfMeasurePoints; j++) {
+        outFile << MP_name[j] << " " << MP_k[j] << " " << MP_level[j] << std::endl;
+        // std::cout << MP_name[j] << "\t" << MP_k[j] << "\t" << MP_level[j] << std::endl;
+    }
+    // std::cout << "Done writing new file..." << std::endl;
+    outFile.close();
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+}
+
 std::string chooseVariation()
 {
     switch (variant) {
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/BoundaryQs.cpp b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/BoundaryQs.cpp
index a0539073347f999e680a9d96ac3d02fee65d7bec..ab266e499fbe9aae16eef29c8bd35b67019323f3 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/BoundaryQs.cpp
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/BoundaryQs.cpp
@@ -166,13 +166,26 @@ unsigned int BoundaryQs::getLevel()
 }
 
 
-void BoundaryQs::setValues(real** q27, unsigned int level) const
+void BoundaryQs::setValuesInVector(std::vector<std::vector<std::vector<real>>> &q27, unsigned int level) const 
+{
+    for (std::size_t column = 0; column < values[level].size(); column++)
+        for (std::size_t index = 0; index < values[level][column].size(); index++)
+            q27[level][column].push_back(values[level][column][index]);
+}
+
+void BoundaryQs::setValues(real **q27, unsigned int level) const
 {
 	for (std::size_t column = 0; column < values[level].size(); column++)
 		for (std::size_t index = 0; index < values[level][column].size(); index++)
 			q27[column][index] = values[level][column][index];
 }
 
+void BoundaryQs::setIndexInVector(std::vector<std::vector<int>> &data, unsigned int level) const 
+{
+    for (std::size_t index = 0; index < indices[level].size(); index++)
+        data[level].push_back(indices[level][index]);
+}
+
 void BoundaryQs::setIndex(int *data, unsigned int level) const
 {
 	for (std::size_t index = 0; index < indices[level].size(); index++)
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/BoundaryQs.h b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/BoundaryQs.h
index daefd2cac6bce4b35e52d915b265be5fd6bd4f99..804051824b5c8ca01809c78e58034eb06d45874e 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/BoundaryQs.h
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/BoundaryQs.h
@@ -32,7 +32,9 @@ private:
 	void init_Binary();
 
 public:
-	void setIndex(int *indices, unsigned int level) const;
+    void setIndexInVector(std::vector<std::vector<int>> &data, unsigned int level) const;
+    void setValuesInVector(std::vector<std::vector<std::vector<real>>> &q27, unsigned int level) const;
+    void setIndex(int *indices, unsigned int level) const;
 	void setValues(real** q27, unsigned int level) const;
 	void getQs(std::vector<std::vector<std::vector<real> > > &qs);
 	void getIndices(std::vector<std::vector<uint> > &indices);
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.cpp b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.cpp
index 49de1d53a3bc0e64907906e2c87409c34cb8ba85..940dbfa617ccd06f5d7b77527cc78b618062240a 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.cpp
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.cpp
@@ -61,10 +61,11 @@ void GridReader::allocArrays_CoordNeighborGeo()
 	CoordNeighborGeoV coordX(para->getcoordX(), binaer, true);
 	CoordNeighborGeoV coordY(para->getcoordY(), binaer, true);
 	CoordNeighborGeoV coordZ(para->getcoordZ(), binaer, true);
-	neighX = std::shared_ptr<CoordNeighborGeoV>(new CoordNeighborGeoV(para->getneighborX(), binaer, false));
-	neighY = std::shared_ptr<CoordNeighborGeoV>(new CoordNeighborGeoV(para->getneighborY(), binaer, false));
-	neighZ = std::shared_ptr<CoordNeighborGeoV>(new CoordNeighborGeoV(para->getneighborZ(), binaer, false));
-	CoordNeighborGeoV geoV(para->getgeoVec(), binaer, false);
+	neighX   = std::shared_ptr<CoordNeighborGeoV>(new CoordNeighborGeoV(para->getneighborX(),   binaer, false));
+	neighY   = std::shared_ptr<CoordNeighborGeoV>(new CoordNeighborGeoV(para->getneighborY(),   binaer, false));
+	neighZ   = std::shared_ptr<CoordNeighborGeoV>(new CoordNeighborGeoV(para->getneighborZ(),   binaer, false));
+    neighWSB = std::shared_ptr<CoordNeighborGeoV>(new CoordNeighborGeoV(para->getneighborWSB(), binaer, false));
+    CoordNeighborGeoV geoV(para->getgeoVec(), binaer, false);
 
 	uint maxLevel = coordX.getLevel();
 	std::cout << "Number of Level: " << maxLevel + 1 << std::endl;
@@ -81,19 +82,20 @@ void GridReader::allocArrays_CoordNeighborGeo()
 
         cudaMemoryManager->cudaAllocCoord(level);
 		cudaMemoryManager->cudaAllocSP(level);
-        cudaMemoryManager->cudaAllocF3SP(level);
+        //cudaMemoryManager->cudaAllocF3SP(level);
         cudaMemoryManager->cudaAllocNeighborWSB(level);
 
         if (para->getUseWale())
 			cudaMemoryManager->cudaAllocTurbulentViscosity(level);
 
-		coordX.initalCoords(para->getParH(level)->coordX_SP, level);
-		coordY.initalCoords(para->getParH(level)->coordY_SP, level);
-		coordZ.initalCoords(para->getParH(level)->coordZ_SP, level);
-		neighX->initalNeighbors(para->getParH(level)->neighborX_SP, level);
-		neighY->initalNeighbors(para->getParH(level)->neighborY_SP, level);
-		neighZ->initalNeighbors(para->getParH(level)->neighborZ_SP, level);
-		geoV.initalNeighbors(para->getParH(level)->geoSP, level);
+		coordX.initalCoords(      para->getParH(level)->coordX_SP,      level);
+		coordY.initalCoords(      para->getParH(level)->coordY_SP,      level);
+		coordZ.initalCoords(      para->getParH(level)->coordZ_SP,      level);
+		neighX->initalNeighbors(  para->getParH(level)->neighborX_SP,   level);
+		neighY->initalNeighbors(  para->getParH(level)->neighborY_SP,   level);
+		neighZ->initalNeighbors(  para->getParH(level)->neighborZ_SP,   level);
+        neighWSB->initalNeighbors(para->getParH(level)->neighborWSB_SP, level);
+        geoV.initalNeighbors(     para->getParH(level)->geoSP,          level);
         rearrangeGeometry(para.get(), level);
 		setInitalNodeValues(numberOfNodesPerLevel, level);
 
@@ -113,13 +115,26 @@ void GridReader::allocArrays_BoundaryValues()
 	this->setChannelBoundaryCondition();
 	int level = BC_Values[0]->getLevel();
 
+	for (int i = 0; i <= level; i++) {
+        velocityX_BCvalues.push_back(std::vector<real>());
+        velocityY_BCvalues.push_back(std::vector<real>());
+        velocityZ_BCvalues.push_back(std::vector<real>());
+        velocityQs.push_back(std::vector<std::vector<real>>());
+        velocityIndex.push_back(std::vector<int>());
+        for (int j = 0; j < para->getD3Qxx(); j++) {
+            velocityQs[i].push_back(std::vector<real>());
+        }
+    }
+
     for (uint i = 0; i < channelBoundaryConditions.size(); i++)
     {
-        setVelocityValues(i);
-        setPressureValues(i);
-        setOutflowValues(i);
+        if (     this->channelBoundaryConditions[i] == "velocity") { fillVelocityVectors(i); } 
+		else if (this->channelBoundaryConditions[i] == "pressure") { setPressureValues(i); } 
+		else if (this->channelBoundaryConditions[i] == "outflow")  { setOutflowValues(i);  }
     }
 
+	setVelocityValues();
+
 	initalValuesDomainDecompostion(level);
 }
 
@@ -229,37 +244,59 @@ void GridReader::setPressRhoBC(int sizePerLevel, int level, int channelSide) con
 }
 
 
-void GridReader::setVelocityValues(int channelSide) const
+void GridReader::fillVelocityVectors(int channelSide)
 {
-	for (unsigned int level = 0; level <= BC_Values[channelSide]->getLevel(); level++)
+    for (unsigned int level = 0; level <= BC_Values[channelSide]->getLevel(); level++)
 	{
-		int sizePerLevel = BC_Values[channelSide]->getSize(level);
-        setVelocitySizePerLevel(level, sizePerLevel);
+		const int sizePerLevel = BC_Values[channelSide]->getSize(level);
 
 		if (sizePerLevel > 1)
 		{
-			std::cout << "size velocity level " << level << " : " << sizePerLevel << std::endl;
+            // set local vectors per side and level
+            real *veloX_ValuesPerSide = new real[sizePerLevel];
+            real *veloY_ValuesPerSide = new real[sizePerLevel];
+            real *veloZ_ValuesPerSide = new real[sizePerLevel];
+
+            std::cout << "size velocity level " << level << " : " << sizePerLevel << std::endl;
+            BC_Values[channelSide]->setVelocityValues(veloX_ValuesPerSide, veloY_ValuesPerSide, veloZ_ValuesPerSide, level);
+
+            for (int i = 0; i < sizePerLevel; i++) {
+                this->velocityX_BCvalues[level].push_back(veloX_ValuesPerSide[i]);
+                this->velocityY_BCvalues[level].push_back(veloY_ValuesPerSide[i]);
+                this->velocityZ_BCvalues[level].push_back(veloZ_ValuesPerSide[i]);
+            }
+
+			delete[] veloX_ValuesPerSide;
+            delete[] veloY_ValuesPerSide;
+            delete[] veloZ_ValuesPerSide;
+        }        
+	}
 
-            cudaMemoryManager->cudaAllocVeloBC(level);
 
-			setVelocity(level, sizePerLevel, channelSide);
-            cudaMemoryManager->cudaCopyVeloBC(level);
-		}
-	}
 }
 
-void GridReader::setVelocity(int level, int sizePerLevel, int channelSide) const
-{
-	BC_Values[channelSide]->setVelocityValues(para->getParH(level)->Qinflow.Vx, para->getParH(level)->Qinflow.Vy, para->getParH(level)->Qinflow.Vz, level);
+void GridReader::setVelocityValues() { 
+    for (int level = 0; level < (int)(velocityX_BCvalues.size()); level++) {
+        
+		int sizePerLevel = (int) velocityX_BCvalues[level].size();
+        std::cout << "complete size velocity level " << level << " : " << sizePerLevel << std::endl;
+        setVelocitySizePerLevel(level, sizePerLevel);
+        
+		if (sizePerLevel > 1) {
+            cudaMemoryManager->cudaAllocVeloBC(level);
+            setVelocity(level, sizePerLevel);
+			cudaMemoryManager->cudaCopyVeloBC(level);
+        }
+    }
+}
 
+void GridReader::setVelocity(int level, int sizePerLevel) const
+{
 	for (int index = 0; index < sizePerLevel; index++)
 	{
-		para->getParH(level)->Qinflow.Vx[index] = para->getParH(level)->Qinflow.Vx[index] / para->getVelocityRatio();
-		para->getParH(level)->Qinflow.Vy[index] = para->getParH(level)->Qinflow.Vy[index] / para->getVelocityRatio();
-		para->getParH(level)->Qinflow.Vz[index] = para->getParH(level)->Qinflow.Vz[index] / para->getVelocityRatio();
-		//para->getParH(level)->Qinflow.Vx[index] = para->getVelocity();//0.035;
-		//para->getParH(level)->Qinflow.Vy[index] = 0.0;//para->getVelocity();//0.0;
-		//para->getParH(level)->Qinflow.Vz[index] = 0.0;
+        para->getParH(level)->Qinflow.Vx[index] = this->velocityX_BCvalues[level][index] / para->getVelocityRatio();
+        para->getParH(level)->Qinflow.Vy[index] = this->velocityY_BCvalues[level][index] / para->getVelocityRatio();
+        para->getParH(level)->Qinflow.Vz[index] = this->velocityZ_BCvalues[level][index] / para->getVelocityRatio();
 	}
 }
 
@@ -514,12 +551,19 @@ void GridReader::allocArrays_BoundaryQs()
 
 	for (std::size_t i = 0; i < channelBoundaryConditions.size(); i++)
 	{
-		if (this->channelBoundaryConditions[i] == "noSlip") { setNoSlipQs(BC_Qs[i]); }
+		if (     this->channelBoundaryConditions[i] == "noSlip"  ) { setNoSlipQs(BC_Qs[i]);   }
 		else if (this->channelBoundaryConditions[i] == "velocity") { setVelocityQs(BC_Qs[i]); }
-		else if (this->channelBoundaryConditions[i] == "pressure") { setPressQs(BC_Qs[i]); }
-		else if (this->channelBoundaryConditions[i] == "outflow") { setOutflowQs(BC_Qs[i]); }
+		else if (this->channelBoundaryConditions[i] == "pressure") { setPressQs(BC_Qs[i]);    }
+		else if (this->channelBoundaryConditions[i] == "outflow" ) { setOutflowQs(BC_Qs[i]);  }
 	}
 
+	for (int lev = 0; lev < (int)(velocityIndex.size()); lev++) {
+        if (velocityIndex[lev].size() > 1) {
+            copyVectorsToQStruct(velocityQs[lev], velocityIndex[lev], para->getParH(lev)->Qinflow);
+            cudaMemoryManager->cudaCopyVeloBC(lev);
+        }
+    }
+
 	std::shared_ptr<BoundaryQs> obj_geomQ = std::shared_ptr<BoundaryQs>(new BoundaryQs(para->getgeomBoundaryBcQs(), para, "geo", false));
 	if (para->getIsGeo())
 		setGeoQs(obj_geomQ);
@@ -544,15 +588,14 @@ void GridReader::setPressQs(std::shared_ptr<BoundaryQs> boundaryQ) const
 	}
 }
 
-void GridReader::setVelocityQs(std::shared_ptr<BoundaryQs> boundaryQ) const
+void GridReader::setVelocityQs(std::shared_ptr<BoundaryQs> boundaryQ)
 {
 	for (unsigned int level = 0; level <= boundaryQ->getLevel(); level++)
 	{
 		if (hasQs(boundaryQ, level))
 		{
 			this->printQSize("velocity", boundaryQ, level);
-			this->initalQStruct(para->getParH(level)->Qinflow, boundaryQ, level);
-            cudaMemoryManager->cudaCopyVeloBC(level);
+            this->initalVectorForQStruct(velocityQs, velocityIndex, boundaryQ, level);
 		}
 	}
 }
@@ -613,7 +656,34 @@ void GridReader::modifyQElement(std::shared_ptr<BoundaryQs> boundaryQ, unsigned
 /*------------------------------------------------------------------------------------------------*/
 /*---------------------------------------private q methods----------------------------------------*/
 /*------------------------------------------------------------------------------------------------*/
-void GridReader::initalQStruct(QforBoundaryConditions& Q, std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const
+void GridReader::initalVectorForQStruct(std::vector<std::vector<std::vector<real>>> &Qs, std::vector<std::vector<int>> &index, 
+										std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const
+{
+    boundaryQ->setValuesInVector(Qs, level);
+    boundaryQ->setIndexInVector(index, level);
+}
+
+void GridReader::copyVectorsToQStruct(std::vector<std::vector<real>> &Qs,
+                                      std::vector<int> &index, QforBoundaryConditions &Q) const
+{
+    QforBoundaryConditions qTemp;
+    this->setQ27Size(qTemp, Q.q27[0], Q.kQ);
+
+	uint sizeOfValues = (uint)index.size();
+
+	for (int direction = 0; direction < para->getD3Qxx(); direction++) {
+        for (size_t indexQ = 0; indexQ < sizeOfValues; indexQ++) {
+            qTemp.q27[direction][indexQ] = Qs[direction][indexQ]; 
+        }
+    }
+
+    for (size_t indexQ = 0; indexQ < sizeOfValues; indexQ++) {
+        Q.k[indexQ] = index[indexQ];
+    }
+}
+
+void GridReader::initalQStruct(QforBoundaryConditions &Q, std::shared_ptr<BoundaryQs> boundaryQ,
+                               unsigned int level) const
 {
 	QforBoundaryConditions qTemp;
 	this->setQ27Size(qTemp, Q.q27[0], Q.kQ);
@@ -628,7 +698,25 @@ bool GridReader::hasQs(std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level
 
 void GridReader::initalGridInformations()
 {
+    int maxLevel = para->getMaxLevel();
+    std::vector<int> gridX, gridY, gridZ;
+    std::vector<int> distX, distY, distZ;
+
+	for (int i = 0; i <= maxLevel; i++) {
+        gridX.push_back(0);
+        gridY.push_back(0);
+        gridZ.push_back(0);
+        distX.push_back(0);
+        distY.push_back(0);
+        distZ.push_back(0);
+    }
 
+    para->setGridX(gridX);
+    para->setGridY(gridY);
+    para->setGridZ(gridZ);
+    para->setDistX(distX);
+    para->setDistY(distY);
+    para->setDistZ(distZ);
 }
 
 void GridReader::setQ27Size(QforBoundaryConditions &Q, real* QQ, unsigned int sizeQ) const
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.h b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.h
index aafe0d1cb732af656fe861e7cab93af4e0ae078b..f7a4c43062da79d39c43e6822688c51ad55e7442 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.h
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.h
@@ -26,6 +26,13 @@ private:
 	std::shared_ptr<CoordNeighborGeoV> neighX, neighY, neighZ, neighWSB;
 	std::vector<std::shared_ptr<BoundaryValues> > BC_Values;
 
+    std::vector<std::vector<real>> velocityX_BCvalues, velocityY_BCvalues, velocityZ_BCvalues;
+    std::vector<std::vector<std::vector<real>>> velocityQs;
+    std::vector<std::vector<int>> velocityIndex;
+
+    std::vector<std::vector<real>> pressureBCvalues;
+    std::vector<std::vector<real>> outflowBCvalues;
+
 public:
 	GridReader(FILEFORMAT format, std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemoryManager> cudaManager);
     ~GridReader();
@@ -50,21 +57,27 @@ 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 fillVelocityVectors(int channelSide);
+    void setVelocityValues();
+	void setVelocity(int level, int sizePerLevel) const;
 
 	void setOutflowValues(int channelSide) const;
 	void setOutflow(int level, int sizePerLevel, int channelSide) const;
 
 
-	void setPressQs(std::shared_ptr<BoundaryQs> boundaryQ) const;
-	void setVelocityQs(std::shared_ptr<BoundaryQs> boundaryQ) const;
+	//void fillVelocityQVectors(int channelSide);
+    void setPressQs(std::shared_ptr<BoundaryQs> boundaryQ) const;
+	void setVelocityQs(std::shared_ptr<BoundaryQs> boundaryQ);
 	void setOutflowQs(std::shared_ptr<BoundaryQs> boundaryQ) const;
 	void setNoSlipQs(std::shared_ptr<BoundaryQs> boundaryQ) const;
 	void setGeoQs(std::shared_ptr<BoundaryQs> boundaryQ) const;
 	void modifyQElement(std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const;
 
-	void initalQStruct(QforBoundaryConditions& Q, std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const;
+	void initalVectorForQStruct(std::vector<std::vector<std::vector<real>>> &Qs, std::vector<std::vector<int>> &index,
+                                std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const;
+    void copyVectorsToQStruct(std::vector<std::vector<real>> &Qs, std::vector<int> &index,
+                              QforBoundaryConditions &Q) const;
+    void initalQStruct(QforBoundaryConditions &Q, std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const;
 	void printQSize(std::string bc, std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const;
 	void setSizeNoSlip(std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const;
 	void setSizeGeoQs(std::shared_ptr<BoundaryQs> boundaryQ, unsigned int level) const;
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index 42b39be681fffac655e9ee8fff8d04cda79f5783..f7320f4d5e3dd22ea59ab9f03fd1d8547597e5ee 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -81,6 +81,8 @@ void Simulation::init(SPtr<Parameter> para, SPtr<GridProvider> gridProvider, std
    this->para = para;
 
    devCheck(comm->mapCudaDevice(para->getMyID(), para->getNumprocs(), para->getDevices(), para->getMaxDev()));
+   
+   para->initLBMSimulationParameter();
 
    gridProvider->allocAndCopyForcing();
    gridProvider->allocAndCopyQuadricLimiters();
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
index ffabc467d049bf9b93e6b199dc92ad90242031ff..f394ea968a7362264fa77071d8c90b26b4d01348 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -50,7 +50,7 @@ Parameter::Parameter(const vf::basics::ConfigurationFile &configData, int number
     ic.myid = myId;
 
     readConfigData(configData);
-    initLBMSimulationParameter();
+    //initLBMSimulationParameter();
 }
 
 void Parameter::readConfigData(const vf::basics::ConfigurationFile &configData)
@@ -373,12 +373,7 @@ void Parameter::readConfigData(const vf::basics::ConfigurationFile &configData)
         this->setDoRestart(configData.getValue<bool>("DoRestart"));
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     if (configData.contains("NOGL"))
-    {
-        maxlevel = configData.getValue<int>("NOGL") - 1;
-        fine = maxlevel;
-    }
-    parH.resize(maxlevel + 1);
-    parD.resize(maxlevel + 1);
+        setMaxLevel(configData.getValue<int>("NOGL"));
 
     this->setGridX(std::vector<int>(this->getMaxLevel() + 1, 32));
     this->setGridY(std::vector<int>(this->getMaxLevel() + 1, 32));
@@ -622,6 +617,9 @@ void Parameter::setD3Qxx(int d3qxx)
 void Parameter::setMaxLevel(int maxlevel)
 {
     this->maxlevel = maxlevel-1;
+    this->fine     = this->maxlevel;
+    parH.resize(this->maxlevel + 1);
+    parD.resize(this->maxlevel + 1);
 }
 void Parameter::setParticleBasicLevel(int pbl)
 {
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
index c95a79771a4433036a695b3990d306fca6709a3c..60beb9f637644d0654732b1f2ece5df43ae1f456 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
@@ -314,6 +314,7 @@ class VIRTUALFLUIDS_GPU_EXPORT Parameter
 {
 public:
     Parameter(const vf::basics::ConfigurationFile &configData, int numberOfProcesses, int myId);
+    void initLBMSimulationParameter();
 
     std::shared_ptr<LBMSimulationParameter> getParH(int level);
     std::shared_ptr<LBMSimulationParameter> getParD(int level);
@@ -750,11 +751,10 @@ public:
     void setInitialCondition(std::function<void(real, real, real, real &, real &, real &, real &)> initialCondition);
     std::function<void(real, real, real, real &, real &, real &, real &)> &getInitialCondition();
 
-    std::vector<std::shared_ptr<LBMSimulationParameter>> parH;
-    std::vector<std::shared_ptr<LBMSimulationParameter>> parD;
+    std::vector<std::shared_ptr<LBMSimulationParameter>> parH = std::vector<std::shared_ptr<LBMSimulationParameter>>(1);
+    std::vector<std::shared_ptr<LBMSimulationParameter>> parD = std::vector<std::shared_ptr<LBMSimulationParameter>>(1);
 private:
     void readConfigData(const vf::basics::ConfigurationFile &configData);
-    void initLBMSimulationParameter();
 
     bool compOn { false };
     bool diffOn { false };