diff --git a/CMake/cmake_config_files/MULE.config.cmake b/CMake/cmake_config_files/MULE.config.cmake
index e4ccd6f5785dc10b52edf0ad860eac3ba43dcfbd..91e788d2c1f0fc5ad4176408fe2aa54d9aca799d 100644
--- a/CMake/cmake_config_files/MULE.config.cmake
+++ b/CMake/cmake_config_files/MULE.config.cmake
@@ -2,4 +2,4 @@ SET(CMAKE_CUDA_ARCHITECTURES "75")
 
 list(APPEND USER_APPS "apps/gpu/LBM/ActuatorLine")
 list(APPEND USER_APPS "apps/gpu/LBM/BoundaryLayer")
-list(APPEND USER_APPS "apps/gpu/LBM/SphereScaling")
\ No newline at end of file
+list(APPEND USER_APPS "apps/gpu/LBM/SphereScaling")
diff --git a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
index 7dc4d5096d7e3d63197a488c1c936dabb61e85ae..8dc9ea853a7658887959bedf6ad564ddf43682a9 100644
--- a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
+++ b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
@@ -248,7 +248,7 @@ void multipleLevel(const std::string& configPath)
 
     gridBuilder->addCoarseGrid( xGridMin,  0.0,  0.0,
                                 xGridMax,  L_y,  L_z, dx);
-    if(false)// Add refinement
+    if(true)// Add refinement
     {
         gridBuilder->setNumberOfLayers(4,0);
         real xMaxRefinement = readPrecursor? xGridMax-H: xGridMax;   //Stop refinement some distance before outlet if domain ist not periodic
diff --git a/pythonbindings/src/gpu/submodules/grid_generator.cpp b/pythonbindings/src/gpu/submodules/grid_generator.cpp
index d46938d081006daabec9bfa27cdc36db017a6bf9..a10b48cc7b7048c1522a783542dadf66010c84b3 100644
--- a/pythonbindings/src/gpu/submodules/grid_generator.cpp
+++ b/pythonbindings/src/gpu/submodules/grid_generator.cpp
@@ -1,4 +1,5 @@
 #include <pybind11/pybind11.h>
+#include "gpu/GridGenerator/utilities/communication.h"
 #include "gpu/GridGenerator/geometries/Object.h"
 #include "gpu/GridGenerator/geometries/BoundingBox/BoundingBox.h"
 #include "gpu/GridGenerator/geometries/Conglomerate/Conglomerate.h"
@@ -17,10 +18,19 @@ namespace grid_generator
     {  
         py::module gridGeneratorModule = parentModule.def_submodule("grid_generator");
 
+        //TODO:
+        // py::enum_<CommunicationDirections>(gridGeneratorModule, "CommunicationDirections")
+        // .value("MX", CommunicationDirections::MX)
+        // .value("PX", CommunicationDirections::PX)
+        // .value("MY", CommunicationDirections::MY)
+        // .value("PY", CommunicationDirections::PY)
+        // .value("MZ", CommunicationDirections::MZ)
+        // .value("PZ", CommunicationDirections::PZ);
+
         py::class_<GridFactory, std::shared_ptr<GridFactory>>(gridGeneratorModule, "GridFactory")
         .def_static("make", &GridFactory::make, py::return_value_policy::reference);
 
-        py::class_<BoundingBox>(gridGeneratorModule, "BoundingBox")
+        py::class_<BoundingBox, std::shared_ptr<BoundingBox>>(gridGeneratorModule, "BoundingBox")
         .def(py::init<real, real, real, real, real, real>(), py::arg("min_x"), py::arg("max_x"), py::arg("min_y"), py::arg("max_y"), py::arg("min_z"), py::arg("max_z"));
 
         py::class_<Object, std::shared_ptr<Object>>(gridGeneratorModule, "Object");
@@ -60,7 +70,11 @@ namespace grid_generator
         .def("add_geometry", py::overload_cast<Object*>(&MultipleGridBuilder::addGeometry), py::arg("solid_object"))
         .def("add_geometry", py::overload_cast<Object*, uint>(&MultipleGridBuilder::addGeometry), py::arg("solid_object"), py::arg("level"))
         .def("get_number_of_levels", &MultipleGridBuilder::getNumberOfLevels)
-        .def("build_grids", &MultipleGridBuilder::buildGrids, py::arg("lbm_or_gks"), py::arg("enable_thin_walls"));
+        .def("build_grids", &MultipleGridBuilder::buildGrids, py::arg("lbm_or_gks"), py::arg("enable_thin_walls"))
+        .def("set_subdomain_box", &MultipleGridBuilder::setSubDomainBox, py::arg("bounding_box"))
+        .def("find_communication_indices", &MultipleGridBuilder::findCommunicationIndices)
+        .def("set_communication_process", &MultipleGridBuilder::setCommunicationProcess)
+        .def("set_number_of_layers", &MultipleGridBuilder::setNumberOfLayers, py::arg("number_of_layers_fine"), py::arg("number_of_layers_between_levels"));
 
         return gridGeneratorModule;
     }
diff --git a/src/gpu/GridGenerator/VelocitySetter/VelocitySetter.cpp b/src/gpu/GridGenerator/VelocitySetter/VelocitySetter.cpp
index 980475d8fc757252b98255e1d1bae0bf2bdd720c..85df8ecc8fbb97c53c3554399c2802867cf1b20b 100644
--- a/src/gpu/GridGenerator/VelocitySetter/VelocitySetter.cpp
+++ b/src/gpu/GridGenerator/VelocitySetter/VelocitySetter.cpp
@@ -366,7 +366,10 @@ void VTKReader::fillArrays(std::vector<real>& coordsY, std::vector<real>& coords
         }
 
         if(!foundAll)
+        {
+            VF_LOG_CRITICAL("Found no matching precursor neighbors for grid point at y={}, z={} \n", posY, posZ);
             throw std::runtime_error("VTKReader::fillArrays(): Did not find neighbors in the VelocityFileCollection for all points");
+        }
     }
 
     if(perfect_match)
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu
index ef74653d6955c0613f0c0388039c7ecce742dcf5..3e553ad68f4066d613037438c94e2a640c1a78fc 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.cu
@@ -12,7 +12,8 @@
 #include "Parameter/CudaStreamManager.h"
 #include "DataStructureInitializer/GridProvider.h"
 #include "GPU/CudaMemoryManager.h"
-#include "lbm/constants/NumericConstants.h"
+#include <lbm/constants/NumericConstants.h>
+#include <logger/Logger.h>
 
 using namespace vf::lbm::constant;
 
@@ -273,11 +274,10 @@ void ActuatorFarm::init(Parameter* para, GridProvider* gridProvider, CudaMemoryM
     this->initBladeIndices(para, cudaMemoryManager);
     this->initBladeVelocities(cudaMemoryManager);
     this->initBladeForces(cudaMemoryManager);    
-    this->initBoundingSphere(para, cudaMemoryManager);  
+    this->initBoundingSpheres(para, cudaMemoryManager);  
     this->streamIndex = 0;
 }
 
-
 void ActuatorFarm::interact(Parameter* para, CudaMemoryManager* cudaMemoryManager, int level, unsigned int t)
 {
     if (level != this->level) return;
@@ -479,9 +479,9 @@ void ActuatorFarm::initBladeIndices(Parameter* para, CudaMemoryManager* cudaMemo
     cudaMemoryManager->cudaCopyBladeIndicesHtoD(this);
 }
 
-void ActuatorFarm::initBoundingSphere(Parameter* para, CudaMemoryManager* cudaMemoryManager)
+void ActuatorFarm::initBoundingSpheres(Parameter* para, CudaMemoryManager* cudaMemoryManager)
 {
-    std::vector<int> nodesInSphere;
+    std::vector<int> nodesInSpheres;
 
     for(uint turbine=0; turbine<this->numberOfTurbines; turbine++)
     {
@@ -490,22 +490,35 @@ void ActuatorFarm::initBoundingSphere(Parameter* para, CudaMemoryManager* cudaMe
         real posX = this->turbinePosXH[turbine];
         real posY = this->turbinePosYH[turbine];
         real posZ = this->turbinePosZH[turbine];
-            
+
         real sphereRadiusSqrd = sphereRadius*sphereRadius;
+            
+        uint minimumNumberOfNodesPerSphere = (uint)(c4o3*cPi*pow(sphereRadius-this->deltaX, c3o1)/pow(this->deltaX, c3o1));
+        uint nodesInThisSphere = 0;
 
         for (uint j = 1; j <= para->getParH(this->level)->numberOfNodes; j++)
         {
             const real distX = para->getParH(this->level)->coordinateX[j]-posX;
             const real distY = para->getParH(this->level)->coordinateY[j]-posY;
             const real distZ = para->getParH(this->level)->coordinateZ[j]-posZ;
-            if(distSqrd(distX,distY,distZ) < sphereRadiusSqrd) nodesInSphere.push_back(j);
+            if(distSqrd(distX,distY,distZ) < sphereRadiusSqrd) 
+            {
+                nodesInSpheres.push_back(j);
+                nodesInThisSphere++;
+            }
+        }
+
+        if(nodesInThisSphere<minimumNumberOfNodesPerSphere)
+        {
+            VF_LOG_CRITICAL("Found only {} nodes in bounding sphere of turbine no. {}, expected at least {}!", nodesInThisSphere, turbine, minimumNumberOfNodesPerSphere);
+            throw std::runtime_error("ActuatorFarm::initBoundingSpheres: Turbine bounding sphere partially out of domain.");
         }
     }
 
-    this->numberOfIndices = uint(nodesInSphere.size());
+    this->numberOfIndices = uint(nodesInSpheres.size());
 
     cudaMemoryManager->cudaAllocSphereIndices(this);
-    std::copy(nodesInSphere.begin(), nodesInSphere.end(), this->boundingSphereIndicesH);
+    std::copy(nodesInSpheres.begin(), nodesInSpheres.end(), this->boundingSphereIndicesH);
     cudaMemoryManager->cudaCopySphereIndicesHtoD(this);
 }
 
diff --git a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.h b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.h
index d5b089177e1c994036cbe95c7ec6b5e08b6abb15..c1dd1add5f508cd1dac29cd9a8672ceb47e86c6d 100644
--- a/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.h
+++ b/src/gpu/VirtualFluids_GPU/PreCollisionInteractor/ActuatorFarm.h
@@ -4,6 +4,9 @@
 #include "PreCollisionInteractor.h"
 #include "PointerDefinitions.h"
 #include "lbm/constants/NumericConstants.h"
+#include <stdexcept>
+
+using namespace vf::lbm::constant;
 
 class Parameter;
 class GridProvider;
@@ -36,6 +39,9 @@ public:
         this->deltaX = _deltaX*exp2(-this->level);
         this->invEpsilonSqrd = 1/(epsilon*epsilon);
         this->invDeltaX = c1o1/this->deltaX;
+     
+        if(this->epsilon<this->deltaX)
+            throw std::runtime_error("ActuatorFarm::ActuatorFarm: epsilon needs to be larger than dx!");
     }
 
     virtual  ~ActuatorFarm()
@@ -141,7 +147,7 @@ public:
 
 private:
     void initTurbineGeometries(CudaMemoryManager* cudaManager);
-    void initBoundingSphere(Parameter* para, CudaMemoryManager* cudaManager);
+    void initBoundingSpheres(Parameter* para, CudaMemoryManager* cudaManager);
     void initBladeCoords(CudaMemoryManager* cudaManager);
     void initBladeVelocities(CudaMemoryManager* cudaManager);
     void initBladeForces(CudaMemoryManager* cudaManager);
diff --git a/src/lbm/constants/NumericConstants.h b/src/lbm/constants/NumericConstants.h
index 1a1350604bf23936cfe091a0291a0f3392697315..b61e994e726e0b5cd85d8130bd6548ccebdc27f6 100644
--- a/src/lbm/constants/NumericConstants.h
+++ b/src/lbm/constants/NumericConstants.h
@@ -49,6 +49,7 @@ static constexpr double c99o100 = 0.99;
 static constexpr double c1o126 = 0.007936507936508;
 static constexpr double c1o216 = 0.004629629629630;
 static constexpr double c5o4 = 1.25;
+static constexpr double c4o3 = 1.333333333333333;
 static constexpr double c9o4 = 2.25;
 static constexpr double c5o2 = 2.5;
 static constexpr double c9o2 = 4.5;
@@ -164,6 +165,7 @@ static constexpr float c99o100 = 0.99f;
 static constexpr float c1o126 = (1.0f / 126.0f);
 static constexpr float c1o216 = (1.0f / 216.0f);
 static constexpr float c5o4 = 1.25f;
+static constexpr float c4o3 = (4.0f / 3.0f);
 static constexpr float c9o4 = 2.25f;
 static constexpr float c5o2 = 2.5f;
 static constexpr float c9o2 = 4.5f;