diff --git a/apps/gpu/LBM/MusselOyster/MusselOyster.cpp b/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
index 702b8851f1c2db5fd8b68cd041176c228e05811f..2b9ce1f112882daff7a22e55b8206ab0b980c34c 100644
--- a/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
+++ b/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
@@ -222,7 +222,7 @@ void multipleLevel(const std::string& configPath)
 
             gridBuilder->buildGrids(LBM, true); // buildGrids() has to be called before setting the BCs!!!!
 
-            gridBuilder->findGeoFluidNodes();
+            gridBuilder->findFluidNodes(useStreams);
 
             if (generatePart == 0) {
                 gridBuilder->findCommunicationIndices(CommunicationDirections::PY, LBM);
@@ -262,7 +262,7 @@ void multipleLevel(const std::string& configPath)
             gridBuilder->setPeriodicBoundaryCondition(false, false, true);
 
             gridBuilder->buildGrids(LBM, true); // buildGrids() has to be called before setting the BCs!!!!
-            gridBuilder->findGeoFluidNodes();
+            gridBuilder->findFluidNodes(useStreams);
 
             //////////////////////////////////////////////////////////////////////////
             gridBuilder->setVelocityBoundaryCondition(SideType::PY, vxLB, 0.0, 0.0);
diff --git a/src/gpu/GridGenerator/grid/Field.cu b/src/gpu/GridGenerator/grid/Field.cu
index 6272a3a6e324395de0a22705ac6e9a61b2879e44..eff280f5aea2a8fb7a1cc751d24bd13a790d163d 100644
--- a/src/gpu/GridGenerator/grid/Field.cu
+++ b/src/gpu/GridGenerator/grid/Field.cu
@@ -67,6 +67,12 @@ HOSTDEVICE bool Field::isFluid(uint index) const
     return type == FLUID || type == FLUID_CFC || type == FLUID_CFF || type == FLUID_FCC || type == FLUID_FCF || isBoundaryConditionNode(index);
 }
 
+HOSTDEVICE bool Field::isFluidBulk(uint index) const
+{
+    const char type = field[index];
+    return type == FLUID || type == FLUID_CFC || type == FLUID_CFF || type == FLUID_FCC || type == FLUID_FCF;
+}
+
 HOSTDEVICE bool Field::isInvalidSolid(uint index) const
 {
     return field[index] == INVALID_SOLID;
diff --git a/src/gpu/GridGenerator/grid/Field.h b/src/gpu/GridGenerator/grid/Field.h
index f79797afab158f0f80d515f2aa38ed8b91438dc2..92b420d6756bdc71012305b483a9efb4641713fc 100644
--- a/src/gpu/GridGenerator/grid/Field.h
+++ b/src/gpu/GridGenerator/grid/Field.h
@@ -22,6 +22,7 @@ public:
     HOSTDEVICE bool isCoarseToFineNode(uint index) const;
     HOSTDEVICE bool isFineToCoarseNode(uint index) const;
 	HOSTDEVICE bool isFluid(uint index) const;
+    HOSTDEVICE bool isFluidBulk(uint index) const;
 	HOSTDEVICE bool isInvalidSolid(uint index) const;
 	HOSTDEVICE bool isQ(uint index) const;
     HOSTDEVICE bool isBoundaryConditionNode(uint index) const;
diff --git a/src/gpu/GridGenerator/grid/Grid.h b/src/gpu/GridGenerator/grid/Grid.h
index 3a70c16add28cac2f717316ce53b1f735d04d7ec..855196dba6070cd3eeeea282cf13003e9228b89b 100644
--- a/src/gpu/GridGenerator/grid/Grid.h
+++ b/src/gpu/GridGenerator/grid/Grid.h
@@ -142,9 +142,10 @@ public:
     virtual void repairCommunicationInices(int direction) = 0;
 
     // needed for CUDA Streams MultiGPU
-    virtual void findFluidNodeIndices() = 0;
-    virtual uint getNumberOfFluidNodes() const = 0;
+    virtual void findFluidNodeIndices(bool onlyBulk) = 0;
+    virtual uint getNumberOfFluidNodes() const = 0;;
     virtual void getFluidNodeIndices(uint *fluidNodeIndices) const = 0;
+
 };
 
 #endif
diff --git a/src/gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp b/src/gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp
index f6e2de17da0520008dbf3998e57e9acbc75c25ca..21570f15d8adb01f18449f788e7a11288aa9c460 100644
--- a/src/gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp
+++ b/src/gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp
@@ -600,11 +600,11 @@ void MultipleGridBuilder::findCommunicationIndices(int direction, LbmOrGks lbmOr
     *logging::out << logging::Logger::INFO_HIGH << "Done with findCommunicationIndices()\n";
 }
 
-void MultipleGridBuilder::findGeoFluidNodes()
+void MultipleGridBuilder::findFluidNodes(bool onlyBulk)
 {
     *logging::out << logging::Logger::INFO_HIGH << "Start findFluidNodes()\n";
     for (uint i = 0; i < grids.size(); i++)
-        grids[i]->findFluidNodeIndices();
+        grids[i]->findFluidNodeIndices(onlyBulk);
     *logging::out << logging::Logger::INFO_HIGH << "Done with findFluidNodes()\n";
 }
 
diff --git a/src/gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h b/src/gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h
index 8039aedd1619ab59f7a9968007b4ab206d679b7b..ad72e1adb1125f00a4240c3828a869f36059a8ba 100644
--- a/src/gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h
+++ b/src/gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h
@@ -86,7 +86,7 @@ public:
     GRIDGENERATOR_EXPORT void findCommunicationIndices( int direction, LbmOrGks lbmOrGks );
 
     // needed for CUDA Streams MultiGPU
-    void findGeoFluidNodes();
+    void findFluidNodes(bool onlyBulk);
 };
 
 #endif
diff --git a/src/gpu/GridGenerator/grid/GridImp.cu b/src/gpu/GridGenerator/grid/GridImp.cu
index d67d319dc872b4f92e4b6bec63d9e9d52b9913c1..8bc404c3aa220e6aa55d9badd86ce138b1705ab8 100644
--- a/src/gpu/GridGenerator/grid/GridImp.cu
+++ b/src/gpu/GridGenerator/grid/GridImp.cu
@@ -858,20 +858,19 @@ CUDA_HOST void GridImp::updateSparseIndices()
     sparseSize = size - removedNodes;
 }
 
-CUDA_HOST void GridImp::findFluidNodeIndices() 
+CUDA_HOST void GridImp::findFluidNodeIndices(bool onlyBulk) 
 {
     this->fluidNodeIndices.clear();
     for (uint index = 0; index < this->size; index++) {
         int sparseIndex = this->getSparseIndex(index);
         if (sparseIndex == -1)
             continue;
-
-        if (this->field.isFluid(index))
+        if ((onlyBulk && this->field.isFluidBulk(index)) || (!onlyBulk && this->field.isFluid(index)))
             this->fluidNodeIndices.push_back((uint)sparseIndex + 1);   // + 1 for numbering shift between GridGenerator and VF_GPU
     }
-    this->numberOfFluidNodes = (uint)this->fluidNodeIndices.size();
 }
 
+
 HOSTDEVICE void GridImp::setNeighborIndices(uint index)
 {
     real x, y, z;
@@ -1762,7 +1761,7 @@ HOSTDEVICE uint GridImp::getSparseSize() const
 }
 
 HOSTDEVICE uint GridImp::getNumberOfFluidNodes() const { 
-    return this->numberOfFluidNodes; 
+    return (uint)this->fluidNodeIndices.size(); 
 }
 
 HOSTDEVICE Field GridImp::getField() const
@@ -1958,11 +1957,10 @@ CUDA_HOST void GridImp::getNodeValues(real *xCoords, real *yCoords, real *zCoord
 }
 
 CUDA_HOST void GridImp::getFluidNodeIndices(uint *fluidNodeIndices) const { 
-  for (uint nodeNumber = 0; nodeNumber < this->numberOfFluidNodes; nodeNumber++)
+    for (uint nodeNumber = 0; nodeNumber < (uint)this->fluidNodeIndices.size(); nodeNumber++)
         fluidNodeIndices[nodeNumber] = this->fluidNodeIndices[nodeNumber];
 }
 
-
 void GridImp::print() const
 {
     printf("min: (%2.4f, %2.4f, %2.4f), max: (%2.4f, %2.4f, %2.4f), size: %d, delta: %2.4f\n", startX, startY, startZ,
diff --git a/src/gpu/GridGenerator/grid/GridImp.h b/src/gpu/GridGenerator/grid/GridImp.h
index 59569809348cda4e702bc59ae38eca871195c9b8..33953ab2a82f3b47993961941780cfcb0f4e8443 100644
--- a/src/gpu/GridGenerator/grid/GridImp.h
+++ b/src/gpu/GridGenerator/grid/GridImp.h
@@ -87,7 +87,6 @@ private:
     int *sparseIndices;
 
     std::vector<uint> fluidNodeIndices;
-    uint numberOfFluidNodes;
 
 	uint *qIndices;     //maps from matrix index to qIndex
 	real *qValues;
@@ -316,7 +315,7 @@ public:
 
     void repairCommunicationInices(int direction) override;
 
-    void findFluidNodeIndices() override;
+    void findFluidNodeIndices(bool onlyBulk) override;
     uint getNumberOfFluidNodes() const override;
     CUDA_HOST void getFluidNodeIndices(uint *fluidNodeIndices) const override;