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 c322f48b98bb5a8107e83802d161725b71e888bf..15250b31e25e1f45a09d976eb6e721458bcba6d2 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 (size_t 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 (size_t 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 (size_t 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 (size_t level = 0; level < velocityX_BCvalues.size(); level++) {
+        
+		int sizePerLevel = 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();
 	}
 }
 
@@ -511,15 +548,23 @@ void GridReader::allocArrays_BoundaryQs()
 
 	std::vector<std::shared_ptr<BoundaryQs> > BC_Qs(channelDirections.size());
 	this->makeReader(BC_Qs, para);
+    int level = BC_Values[0]->getLevel();
 
 	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 (size_t lev = 0; lev < 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 +589,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 +657,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 = index.size();
+
+	for (size_t 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);
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;