diff --git a/apps/gpu/LBM/MusselOyster/MusselOyster.cpp b/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
index 7f0aa584f1928f670f4cfe71ec7e056ac92d2296..702b8851f1c2db5fd8b68cd041176c228e05811f 100644
--- a/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
+++ b/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
@@ -110,6 +110,7 @@ void multipleLevel(const std::string& configPath)
 
     bool useGridGenerator = true;
     bool useMultiGPU = false;
+    bool useStreams= true;
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -142,9 +143,10 @@ void multipleLevel(const std::string& configPath)
     para->setFName(para->getOutputPath() + "/" + para->getOutputPrefix());
     para->setPrintFiles(true);
 
-    para->setMaxLevel(2);
+    para->setMaxLevel(1);
 
     //para->setMainKernel("CumulantK17CompChim");
+    para->useStreams = useStreams;
     para->setMainKernel("CumulantK17CompChimSparse");
    *logging::out << logging::Logger::INFO_HIGH << "Kernel: " << para->getMainKernel() << "\n";
 
@@ -252,8 +254,8 @@ void multipleLevel(const std::string& configPath)
 
             gridBuilder->addCoarseGrid(xGridMin, yGridMin, zGridMin, xGridMax, yGridMax, zGridMax, dxGrid);
 
-            gridBuilder->setNumberOfLayers(6, 8);
-            gridBuilder->addGrid(bivalveRef_1_STL, 1);
+            //gridBuilder->setNumberOfLayers(6, 8);
+            //gridBuilder->addGrid(bivalveRef_1_STL, 1);
 
             gridBuilder->addGeometry(bivalveSTL);
 
@@ -347,21 +349,22 @@ void multipleLevel(const std::string& configPath)
         std::cout << ".....level 0: sparse index geoFluid \t" << sparseIndicesFluid[i] << ",    fluid nodes index  \t"
                   << para->getParH(0)->fluidNodeIndices[i] << std::endl;
 
-    uint *geoSP1      = para->getParH(1)->geoSP;
-    uint numGeoFluid1 = 0;
-    std::vector<int> sparseIndicesFluid1;
-    for (uint i = 0; i < para->getParH(1)->size_Mat_SP; i++) {
-        if (geoSP1[i] == GEO_FLUID) {
-            numGeoFluid1++;
-            sparseIndicesFluid1.push_back(i);
-        }
-    }
-    std::cout << ".....geoFluid level 1 " << numGeoFluid1 << ", num fluid nodes (new kernel)  "
-              << para->getParH(1)->numberOfFluidNodes << std::endl;
-
-    for (uint i = 300000; i < 300003; i++)
-        std::cout << ".....level 1: sparse index geoFluid \t" << sparseIndicesFluid[i] << ",    fluid nodes index  \t"
-                  << para->getParH(0)->fluidNodeIndices[i] << std::endl;
+    //// Level 1
+    //uint *geoSP1      = para->getParH(1)->geoSP;
+    //uint numGeoFluid1 = 0;
+    //std::vector<int> sparseIndicesFluid1;
+    //for (uint i = 0; i < para->getParH(1)->size_Mat_SP; i++) {
+    //    if (geoSP1[i] == GEO_FLUID) {
+    //        numGeoFluid1++;
+    //        sparseIndicesFluid1.push_back(i);
+    //    }
+    //}
+    //std::cout << ".....geoFluid level 1 " << numGeoFluid1 << ", num fluid nodes (new kernel)  "
+    //          << para->getParH(1)->numberOfFluidNodes << std::endl;
+
+    //for (uint i = 300000; i < 300003; i++)
+    //    std::cout << ".....level 1: sparse index geoFluid \t" << sparseIndicesFluid[i] << ",    fluid nodes index  \t"
+    //              << para->getParH(0)->fluidNodeIndices[i] << std::endl;
 
 
     sim.run();
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index cbfc3072ac6b538efd36a91a50d53a5c39eb5aa7..08703cdff85435c9e346927983d74a102df15a15 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -103,6 +103,11 @@ void Simulation::init(SPtr<Parameter> para, SPtr<GridProvider> gridProvider, std
    if(para->getMyID() == 0) output.setConsoleOut(true);
    output.clearLogFile();
    //////////////////////////////////////////////////////////////////////////
+   // CUDA streams
+   if(para->useStreams)
+	   para->launchStreams((uint)1);
+   //////////////////////////////////////////////////////////////////////////
+   // 
    //output << para->getNeedInterface().at(0) << "\n";
    //output << para->getNeedInterface().at(1) << "\n";
    //output << para->getNeedInterface().at(2) << "\n";
@@ -1150,6 +1155,10 @@ void Simulation::definePMarea(std::shared_ptr<PorousMedia> pMedia)
 
 void Simulation::free()
 {
+	// Cuda Streams
+    if (para->useStreams)
+        para->terminateStreams();
+
 	//CudaFreeHostMemory
 	for (int lev = para->getCoarse(); lev <= para->getFine(); lev++)
 	{
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
index bf03b915e927ae484ac925eac968c9f5998f93a9..7bf40e27b6968cf115194fa5dffe4ff2790e7407 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -53,6 +53,24 @@ Parameter::Parameter(const vf::basics::ConfigurationFile &configData, int number
     //initLBMSimulationParameter();
 }
 
+void Parameter::launchStreams(uint numberOfStreams)
+{
+    cudaStreams.resize(numberOfStreams);
+    for (cudaStream_t &stream : cudaStreams) {
+        cudaStreamCreate(&stream);
+    }
+}
+
+void Parameter::terminateStreams() {
+    for (cudaStream_t &stream : cudaStreams) {
+        cudaStreamDestroy(stream);
+    }
+}
+
+cudaStream_t& Parameter::getStream(uint streamIndex) { 
+	return cudaStreams[streamIndex]; 
+}
+
 void Parameter::readConfigData(const vf::basics::ConfigurationFile &configData)
 {
    if (configData.contains("NumberOfDevices"))
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
index 80b2fac298f6594f86fe34a1355c8ad594327908..d44c3abd6f061306315283b5fe573be08e5969d8 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
@@ -43,6 +43,8 @@
 
 #include "VirtualFluids_GPU_export.h"
 
+#include <cuda_runtime.h>
+#include <helper_cuda.h>
 
 struct curandStateXORWOW;
 typedef struct curandStateXORWOW curandState;
@@ -755,6 +757,14 @@ public:
 
     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);
+
+    ////////////////////////////////////////////////////////////////////////////
+    // cuda streams
+    bool useStreams = false;
+    void launchStreams(uint numberOfStreams);
+    void terminateStreams();
+    cudaStream_t& getStream(uint streamIndex);
+
 private:
     void readConfigData(const vf::basics::ConfigurationFile &configData);
 
@@ -834,6 +844,8 @@ private:
     std::vector<std::string> possNeighborFilesRecvX, possNeighborFilesRecvY, possNeighborFilesRecvZ;
     bool isNeigborX, isNeigborY, isNeigborZ;
 
+    std::vector<cudaStream_t> cudaStreams;
+
     ////////////////////////////////////////////////////////////////////////////
     // initial condition
     std::function<void(real, real, real, real &, real &, real &, real &)> initialCondition;