diff --git a/CMake/cmake_config_files/ARAGORN.config.cmake b/CMake/cmake_config_files/ARAGORN.config.cmake
index c8432efe045c386174a9e2a04988ed51ed794bf3..d713f02d971024f29d3fb0fd30cfce7585d9dc55 100644
--- a/CMake/cmake_config_files/ARAGORN.config.cmake
+++ b/CMake/cmake_config_files/ARAGORN.config.cmake
@@ -1,14 +1,19 @@
 #################################################################################
 # VirtualFluids MACHINE FILE
 # Responsible: Anna Wellmann
-# OS:          Ubuntu 20.04 (Docker container)
+# OS:          Windows 11
 #################################################################################
 
 set(CMAKE_CUDA_ARCHITECTURES 86)     # Nvidia GeForce RTX 3060
 
+# numerical tests location of the grids
+# SET(PATH_NUMERICAL_TESTS "E:/temp/numericalTests/")
+# list(APPEND VF_COMPILER_DEFINITION "PATH_NUMERICAL_TESTS=${PATH_NUMERICAL_TESTS}")
+
+# add invidual apps here
 set(GPU_APP "apps/gpu/LBM/")
 list(APPEND USER_APPS 
     "${GPU_APP}DrivenCavityMultiGPU"
     "${GPU_APP}SphereScaling"
     # "${GPU_APP}MusselOyster"
-    )
\ No newline at end of file
+    )
diff --git a/CMake/cmake_config_files/ARAGORNUBUNTU.config.cmake b/CMake/cmake_config_files/ARAGORNUBUNTU.config.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..c8432efe045c386174a9e2a04988ed51ed794bf3
--- /dev/null
+++ b/CMake/cmake_config_files/ARAGORNUBUNTU.config.cmake
@@ -0,0 +1,14 @@
+#################################################################################
+# VirtualFluids MACHINE FILE
+# Responsible: Anna Wellmann
+# OS:          Ubuntu 20.04 (Docker container)
+#################################################################################
+
+set(CMAKE_CUDA_ARCHITECTURES 86)     # Nvidia GeForce RTX 3060
+
+set(GPU_APP "apps/gpu/LBM/")
+list(APPEND USER_APPS 
+    "${GPU_APP}DrivenCavityMultiGPU"
+    "${GPU_APP}SphereScaling"
+    # "${GPU_APP}MusselOyster"
+    )
\ No newline at end of file
diff --git a/apps/gpu/LBM/DrivenCavityMultiGPU/DrivenCavityMultiGPU.cpp b/apps/gpu/LBM/DrivenCavityMultiGPU/DrivenCavityMultiGPU.cpp
index 4bbd7f07990c3093309124fb3cc9e9e3852b6574..cd0c73a3ecf8af1a534a53ba1237aea64bfe4b59 100644
--- a/apps/gpu/LBM/DrivenCavityMultiGPU/DrivenCavityMultiGPU.cpp
+++ b/apps/gpu/LBM/DrivenCavityMultiGPU/DrivenCavityMultiGPU.cpp
@@ -67,20 +67,8 @@
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-//  Tesla 03
-// const std::string outPath("E:/temp/DrivenCavityMultiGPUResults/");
-// const std::string gridPath = "D:/STLs/DrivenCavity";
-// const std::string simulationName("DrivenCavityMultiGPU");
-
-// Phoenix
-// const std::string outPath("/work/y0078217/Results/DrivenCavityMultiGPUResults/");
-// const std::string gridPath = "/work/y0078217/Grids/GridDrivenCavityMultiGPU/";
-// const std::string simulationName("DrivenCavityMultiGPU");
-
-//  Aragorn
 const std::string outPath("output/DrivenCavity_Results/");
 const std::string gridPath = "output/DrivenCavity_Results/grid/";
-const std::string simulationName("DrivenCavity");
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -136,16 +124,6 @@ void multipleLevel(std::filesystem::path& configPath)
     const real vyLB        = velocityLB / (real)sqrt(2.0); // LB units
     const real viscosityLB = nx * velocityLB / Re;         // LB units
 
-    VF_LOG_INFO("LB parameters:");
-    VF_LOG_INFO("velocity LB [dx/dt]              = {}", vxLB);
-    VF_LOG_INFO("viscosity LB [dx/dt]             = {}", viscosityLB);
-    VF_LOG_INFO("dxGrid [-]                       = {}\n", dxGrid);
-
-    para->setVelocityLB(velocityLB);
-    para->setViscosityLB(viscosityLB);
-    para->setVelocityRatio(velocity / velocityLB);
-    para->setDensityRatio((real)1.0); // correct value?
-
     para->setInitialCondition([&](real coordX, real coordY, real coordZ, real &rho, real &vx, real &vy, real &vz) {
         rho = (real)1.0;
         vx  = (real)(coordX * velocityLB);
@@ -153,11 +131,17 @@ void multipleLevel(std::filesystem::path& configPath)
         vz  = (real)(coordZ * velocityLB);
     });
 
+    para->setVelocityLB(velocityLB);
+    para->setViscosityLB(viscosityLB);
+    para->setVelocityRatio(velocity / velocityLB);
+    para->setDensityRatio((real)1.0); // correct value?
+
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     para->setCalcDragLift(false);
     para->setUseWale(false);
 
+    if (para->getOutputPath() == "output/") {para->setOutputPath(outPath);}
     para->setOutputPrefix(simulationName);
 
     para->setPrintFiles(true);
@@ -165,7 +149,17 @@ void multipleLevel(std::filesystem::path& configPath)
 
     // para->setMainKernel("CumulantK17CompChim");
     para->setMainKernel("CumulantK17CompChimStream");
-    *logging::out << logging::Logger::INFO_HIGH << "Kernel: " << para->getMainKernel() << "\n";
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    vf::logging::Logger::changeLogPath(para->getOutputPath());
+    VF_LOG_INFO("LB parameters:");
+    VF_LOG_INFO("velocity LB [dx/dt]              = {}", vxLB);
+    VF_LOG_INFO("viscosity LB [dx/dt]             = {}", viscosityLB);
+    VF_LOG_INFO("dxGrid [-]                       = {}\n", dxGrid);
+
+    VF_LOG_INFO("simulation parameters:");
+    VF_LOG_INFO("mainKernel                       = {}\n", para->getMainKernel());
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -559,7 +553,7 @@ int main(int argc, char *argv[])
             std::filesystem::path configPath = __FILE__;
 
             // the config file's default name can be replaced by passing a command line argument
-            std::string configName("config.txt");
+            std::string configName("configDrivenCavityMultiGPU.txt");
             if (argc == 2) {
                 configName = argv[1];
                 std::cout << "Using configFile command line argument: " << configName << std::endl;
diff --git a/apps/gpu/LBM/MusselOyster/MusselOyster.cpp b/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
index aa782d942f8aea811cc59a1d17df2319ef0e5df5..106c9de23f0458b5ace0fcf988e7d1f53b4e71c9 100644
--- a/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
+++ b/apps/gpu/LBM/MusselOyster/MusselOyster.cpp
@@ -64,18 +64,18 @@
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
+// Relative Paths
+const std::string outPath("./output/MusselOysterResults/");
+const std::string gridPathParent = "./output/MusselOysterResults/grid/";
+const std::string stlPath("./stl/MusselOyster/");
+const std::string simulationName("MusselOyster");
+
 // Tesla 03
 // const std::string outPath("E:/temp/MusselOysterResults/");
 // const std::string gridPathParent = "E:/temp/GridMussel/";
 // const std::string stlPath("C:/Users/Master/Documents/MasterAnna/STL/");
 // const std::string simulationName("MusselOyster");
 
-// Aragorn
-const std::string outPath("./output/MusselOysterResults/");
-const std::string gridPathParent = "./output/MusselOysterResults/grid/";
-const std::string stlPath("./stl/MusselOyster/");
-const std::string simulationName("MusselOyster");
-
 // Phoenix
 // const std::string outPath("/work/y0078217/Results/MusselOysterResults/");
 // const std::string gridPathParent = "/work/y0078217/Grids/GridMusselOyster/";
@@ -126,7 +126,6 @@ void multipleLevel(std::filesystem::path &configPath)
         para->useReducedCommunicationAfterFtoC = false;
     }
 
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     std::string bivalveType = "MUSSEL"; // "MUSSEL" "OYSTER"
     std::string gridPath(
@@ -148,6 +147,24 @@ void multipleLevel(std::filesystem::path &configPath)
     para->setViscosityRatio((real)0.058823529);
     para->setDensityRatio((real)998.0);
 
+    // para->setTimestepOut(1000);
+    // para->setTimestepEnd(10000);
+
+    para->setCalcDragLift(false);
+    para->setUseWale(false);
+
+    para->setOutputPrefix(simulationName);
+    if (para->getOutputPath() == "output/") {para->setOutputPath(outPath);}
+
+    para->setPrintFiles(true);
+    std::cout << "Write result files to " << para->getFName() << std::endl;
+
+    para->setUseStreams(useStreams);
+    // para->setMainKernel("CumulantK17CompChim");
+    para->setMainKernel("CumulantK17CompChimStream");
+    
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    
     VF_LOG_INFO("LB parameters:");
     VF_LOG_INFO("velocity LB [dx/dt]              = {}", vxLB);
     VF_LOG_INFO("viscosity LB [dx/dt]             = {}", viscosityLB);
@@ -165,22 +182,6 @@ void multipleLevel(std::filesystem::path &configPath)
     VF_LOG_INFO("bivalveType                      = {}", bivalveType);
     VF_LOG_INFO("mainKernel                       = {}\n", para->getMainKernel());
 
-    // para->setTimestepOut(1000);
-    // para->setTimestepEnd(10000);
-
-    para->setCalcDragLift(false);
-    para->setUseWale(false);
-
-    para->setOutputPrefix(simulationName);
-
-    para->setPrintFiles(true);
-    std::cout << "Write result files to " << para->getFName() << std::endl;
-
-    para->setUseStreams(useStreams);
-    // para->setMainKernel("CumulantK17CompChim");
-    para->setMainKernel("CumulantK17CompChimStream");
-    *logging::out << logging::Logger::INFO_HIGH << "Kernel: " << para->getMainKernel() << "\n";
-
     //////////////////////////////////////////////////////////////////////////
 
     if (useGridGenerator) {
@@ -597,7 +598,7 @@ int main(int argc, char *argv[])
             std::filesystem::path configPath = __FILE__;
 
             // the config file's default name can be replaced by passing a command line argument
-            std::string configName("config.txt");
+            std::string configName("configMusselOyster.txt");
             if (argc == 2) {
                 configName = argv[1];
                 std::cout << "Using configFile command line argument: " << configName << std::endl;
diff --git a/apps/gpu/LBM/SphereGPU/config.txt b/apps/gpu/LBM/SphereGPU/config.txt
index 7d8b3831490ff828a3d5eaaa047aefeadcd6269e..9d0804ffc0d364e8d75ee16424378f360f07c52d 100644
--- a/apps/gpu/LBM/SphereGPU/config.txt
+++ b/apps/gpu/LBM/SphereGPU/config.txt
@@ -8,7 +8,7 @@
 #informations for Writing
 ##################################################
 Path=output/Sphere/
-Prefix=Sphere
+Prefix=Sphere01
 #WriteGrid=true
 ##################################################
 #informations for reading
diff --git a/apps/gpu/LBM/SphereScaling/SphereScaling.cpp b/apps/gpu/LBM/SphereScaling/SphereScaling.cpp
index 6127491efac1e8f03abf13e06b87e3baaa9d4955..a5f25e3d52dfc03f5247867bb79a6edfbc4a26b4 100644
--- a/apps/gpu/LBM/SphereScaling/SphereScaling.cpp
+++ b/apps/gpu/LBM/SphereScaling/SphereScaling.cpp
@@ -149,12 +149,14 @@ void multipleLevel(std::filesystem::path& configPath)
 
 
     para->setOutputPrefix(simulationName);
+    if (para->getOutputPath() == "output/") {para->setOutputPath(outPath);}
     para->setPrintFiles(true);
-    std::cout << "Write result files to " << para->getFName() << std::endl;
 
     // para->setMainKernel("CumulantK17CompChim");
     para->setMainKernel("CumulantK17CompChimStream");
 
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
     VF_LOG_INFO("LB parameters:");
     VF_LOG_INFO("velocity LB [dx/dt]              = {}", vxLB);
     VF_LOG_INFO("viscosity LB [dx/dt]             = {}", viscosityLB);
diff --git a/src/gpu/GridGenerator/CMakeLists.txt b/src/gpu/GridGenerator/CMakeLists.txt
index 7483d276d08bb66b8c6948689f1dbb13f8846cb6..8102ad3a10b53dded2ba6fe489753f20d1d2ed4f 100644
--- a/src/gpu/GridGenerator/CMakeLists.txt
+++ b/src/gpu/GridGenerator/CMakeLists.txt
@@ -5,4 +5,9 @@ vf_add_tests()
 
 if(NOT MSVC) 
    target_compile_options(GridGenerator PRIVATE "-Wno-strict-aliasing")
+endif()
+
+
+if(BUILD_VF_UNIT_TESTS)
+   target_include_directories(GridGeneratorTests PRIVATE "${VF_ROOT_DIR}/src/basics/")
 endif()
\ No newline at end of file
diff --git a/src/gpu/GridGenerator/grid/Field.cpp b/src/gpu/GridGenerator/grid/Field.cpp
index 1777712107cc630d8b717233c11c224e278f787a..8e68bfdf7e324a13a4c9d31d493580e56a48b6cd 100644
--- a/src/gpu/GridGenerator/grid/Field.cpp
+++ b/src/gpu/GridGenerator/grid/Field.cpp
@@ -38,7 +38,6 @@ using namespace vf::gpu;
 
 Field::Field(uint size) : size(size)
 {
-    
 }
 
 void Field::allocateMemory()
diff --git a/src/gpu/GridGenerator/grid/Field.h b/src/gpu/GridGenerator/grid/Field.h
index db664d8820cce5eaeb01d08ef4d37f1ed7bd7064..bb25b0fc03537b05eaadfa5c3d161f83c1267fae 100644
--- a/src/gpu/GridGenerator/grid/Field.h
+++ b/src/gpu/GridGenerator/grid/Field.h
@@ -71,7 +71,7 @@ public:
     void setFieldEntryToInvalidCoarseUnderFine(uint index);
     void setFieldEntryToInvalidOutOfGrid(uint index);
 
-private:
+protected:
     char *field = nullptr;
     uint size;
 };
diff --git a/src/gpu/GridGenerator/grid/GridImp.cpp b/src/gpu/GridGenerator/grid/GridImp.cpp
index 030f94344817521a6059fcd498cb7fc457a7d862..31bbf3ddc87184846fcb01a3e6631358b6a6f864 100644
--- a/src/gpu/GridGenerator/grid/GridImp.cpp
+++ b/src/gpu/GridGenerator/grid/GridImp.cpp
@@ -993,13 +993,13 @@ void GridImp::setStopperNeighborCoords(uint index)
     real x, y, z;
     this->transIndexToCoords(index, x, y, z);
 
-    if (vf::Math::lessEqual(x + delta, endX) && !this->field.isInvalidOutOfGrid(this->transCoordToIndex(x + delta, y, z)))
+    if (vf::Math::lessEqual(x + delta, endX + (0.5 * delta)) && !this->field.isInvalidOutOfGrid(this->transCoordToIndex(x + delta, y, z)))
         neighborIndexX[index] = getSparseIndex(x + delta, y, z);
 
-    if (vf::Math::lessEqual(y + delta, endY) && !this->field.isInvalidOutOfGrid(this->transCoordToIndex(x, y + delta, z)))
+    if (vf::Math::lessEqual(y + delta, endY + (0.5 * delta)) && !this->field.isInvalidOutOfGrid(this->transCoordToIndex(x, y + delta, z)))
         neighborIndexY[index] = getSparseIndex(x, y + delta, z);
 
-    if (vf::Math::lessEqual(z + delta, endZ) && !this->field.isInvalidOutOfGrid(this->transCoordToIndex(x, y, z + delta)))
+    if (vf::Math::lessEqual(z + delta, endZ + (0.5 * delta)) && !this->field.isInvalidOutOfGrid(this->transCoordToIndex(x, y, z + delta)))
         neighborIndexZ[index] = getSparseIndex(x, y, z + delta);
 
     if (vf::Math::greaterEqual(x - delta, endX) && 
diff --git a/src/gpu/GridGenerator/grid/GridImp.h b/src/gpu/GridGenerator/grid/GridImp.h
index ee30e2b4aaadd737e1fa096eec3b815768ddd0a0..158cba9e67650d65c20e63aa6a35d45f129e2baa 100644
--- a/src/gpu/GridGenerator/grid/GridImp.h
+++ b/src/gpu/GridGenerator/grid/GridImp.h
@@ -109,11 +109,9 @@ private:
     uint sparseSize;
     bool periodicityX = false, periodicityY = false, periodicityZ = false;
 
-    Field field;
     Object* object;
     GridInterface *gridInterface;
 
-    int *neighborIndexX, *neighborIndexY, *neighborIndexZ, *neighborIndexNegative;
     int *sparseIndices;
 
     std::vector<uint> fluidNodeIndices;
@@ -133,6 +131,10 @@ private:
 
     bool enableFixRefinementIntoTheWall;
 
+protected:
+    Field field;
+    int *neighborIndexX, *neighborIndexY, *neighborIndexZ, *neighborIndexNegative;
+
 public:
     void inital(const SPtr<Grid> fineGrid, uint numberOfLayers) override;
     void setOddStart(bool xOddStart, bool yOddStart, bool zOddStart) override;
@@ -155,11 +157,11 @@ public:
     uint transCoordToIndex(const real &x, const real &y, const real &z) const override;
     void transIndexToCoords(uint index, real &x, real &y, real &z) const override;
 
-    virtual void findGridInterface(SPtr<Grid> grid, LbmOrGks lbmOrGks) override;
+    void findGridInterface(SPtr<Grid> grid, LbmOrGks lbmOrGks) override;
 
     void repairGridInterfaceOnMultiGPU(SPtr<Grid> fineGrid) override;
 
-    virtual void limitToSubDomain(SPtr<BoundingBox> subDomainBox, LbmOrGks lbmOrGks) override;
+    void limitToSubDomain(SPtr<BoundingBox> subDomainBox, LbmOrGks lbmOrGks) override;
 
     void freeMemory() override;
 
@@ -277,15 +279,16 @@ public:
     void setNeighborIndices(uint index);
     real getFirstFluidNode(real coords[3], int direction, real startCoord) const override;
     real getLastFluidNode(real coords[3], int direction, real startCoord) const override;
+protected:
+    virtual void setStopperNeighborCoords(uint index);
 private:
-    void setStopperNeighborCoords(uint index);
     void getNeighborCoords(real &neighborX, real &neighborY, real &neighborZ, real x, real y, real z) const;
     real getNeighborCoord(bool periodicity, real endCoord, real coords[3], int direction) const;
     void getNegativeNeighborCoords(real &neighborX, real &neighborY, real &neighborZ, real x, real y, real z) const;
     real getNegativeNeighborCoord(bool periodicity, real endCoord, real coords[3], int direction) const;
     
 
-    int getSparseIndex(const real &expectedX, const real &expectedY, const real &expectedZ) const;
+    virtual int getSparseIndex(const real &expectedX, const real &expectedY, const real &expectedZ) const;
 
     static real getMinimumOnNodes(const real &minExact, const real &decimalStart, const real &delta);
     static real getMaximumOnNodes(const real &maxExact, const real &decimalStart, const real &delta);
diff --git a/src/gpu/GridGenerator/grid/GridImpTest.cpp b/src/gpu/GridGenerator/grid/GridImpTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f5ddb1b01dd88cca7d750017ec328efe02cd92f
--- /dev/null
+++ b/src/gpu/GridGenerator/grid/GridImpTest.cpp
@@ -0,0 +1,258 @@
+#include <array>
+#include <gmock/gmock-matchers.h>
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+#include <memory>
+#include <ostream>
+
+#include "GridImp.h"
+#include "PointerDefinitions.h"
+#include "grid/Field.h"
+#include "grid/GridBuilder/MultipleGridBuilder.h"
+#include "grid/distributions/Distribution.h"
+
+// This test is commented out because it causes a compiler error in Clang 10 --> The bug is fixed in Clang 14 (https://github.com/google/googletest/issues/2271)
+
+// class FieldDouble : public Field
+// {
+// public:
+//     FieldDouble() : Field(1)
+//     {
+//         this->allocateMemory();
+//     };
+
+//     void setToStopper(uint index)
+//     {
+//         this->field[index] = vf::gpu::STOPPER_SOLID;
+//     }
+// };
+
+// class GridImpDouble : public GridImp
+// {
+// public:
+//     std::array<real, 3> coordsOfTestedNode;
+//     GridImpDouble(Object *object, real startX, real startY, real startZ, real endX, real endY, real endZ, real delta,
+//                   Distribution d, uint level)
+//         : GridImp(object, startX, startY, startZ, endX, endY, endZ, delta, d, level)
+//     {
+//         this->neighborIndexX = new int[5];
+//         this->neighborIndexY = new int[5];
+//         this->neighborIndexZ = new int[5];
+//     }
+
+//     static SPtr<GridImpDouble> makeShared(Object *object, real startX, real startY, real startZ, real endX, real endY,
+//                                           real endZ, real delta, Distribution d, uint level)
+//     {
+//         SPtr<GridImpDouble> grid(new GridImpDouble(object, startX, startY, startZ, endX, endY, endZ, delta, d, level));
+//         return grid;
+//     }
+
+//     void transIndexToCoords(uint, real &x, real &y, real &z) const override
+//     {
+//         x = coordsOfTestedNode[0];
+//         y = coordsOfTestedNode[1];
+//         z = coordsOfTestedNode[2];
+//     }
+
+//     uint transCoordToIndex(const real &, const real &, const real &) const override
+//     {
+//         return 0;
+//     }
+
+//     void setStopperNeighborCoords(uint index) override
+//     {
+//         GridImp::setStopperNeighborCoords(index);
+//     }
+
+//     void setField(Field &field)
+//     {
+//         this->field = field;
+//     }
+
+//     MOCK_METHOD(int, getSparseIndex, (const real &x, const real &y, const real &z), (const, override));
+// };
+
+// // This is test is highly dependent on the implementation. Maybe it should be removed :(
+// TEST(GridImp, setStopperNeighborCoords)
+// {
+//     real end = 1.0;
+//     real delta = 0.1;
+
+//     SPtr<GridImpDouble> gridImp =
+//         GridImpDouble::makeShared(nullptr, 0.0, 0.0, 0.0, end, end, end, delta, Distribution(), 0);
+//     FieldDouble field;
+//     field.setToStopper(0);
+//     gridImp->setField(field);
+
+//     gridImp->coordsOfTestedNode = { end - ((real)0.5 * delta), end - ((real)0.5 * delta), end - ((real)0.5 * delta) };
+//     EXPECT_CALL(*gridImp, getSparseIndex).Times(3);
+//     gridImp->setStopperNeighborCoords(0);
+
+//     gridImp->coordsOfTestedNode = { end - ((real)0.51 * delta), end - ((real)0.51 * delta),
+//                                     end - ((real)0.51 * delta) };
+//     EXPECT_CALL(*gridImp, getSparseIndex).Times(3);
+//     gridImp->setStopperNeighborCoords(0);
+//     gridImp->coordsOfTestedNode = { end - ((real)0.99 * delta), end - ((real)0.99 * delta),
+//                                     end - ((real)0.99 * delta) };
+//     EXPECT_CALL(*gridImp, getSparseIndex).Times(3);
+//     gridImp->setStopperNeighborCoords(0);
+
+//     gridImp->coordsOfTestedNode = { end - delta, end - delta, end - delta };
+//     EXPECT_CALL(*gridImp, getSparseIndex).Times(3);
+//     gridImp->setStopperNeighborCoords(0);
+
+//     gridImp->coordsOfTestedNode = { end - ((real)1.01 * delta), end - ((real)1.01 * delta),
+//                                     end - ((real)1.01 * delta) };
+//     EXPECT_CALL(*gridImp, getSparseIndex).Times(3);
+//     gridImp->setStopperNeighborCoords(0);
+
+//     // The grid should not be like this, so this should be fine...
+//     gridImp->coordsOfTestedNode = { end, end, end };
+//     EXPECT_CALL(*gridImp, getSparseIndex).Times(0);
+//     gridImp->setStopperNeighborCoords(0);
+
+//     gridImp->coordsOfTestedNode = { end - ((real)0.25 * delta), end - ((real)0.25 * delta),
+//                                     end - ((real)0.25 * delta) };
+//     EXPECT_CALL(*gridImp, getSparseIndex).Times(0);
+//     gridImp->setStopperNeighborCoords(0);
+// }
+
+std::array<int, 3> countInvalidNeighbors(SPtr<Grid> grid)
+{
+    auto countInvalidX = 0;
+    auto countInvalidY = 0;
+    auto countInvalidZ = 0;
+    for (uint index = 0; index < grid->getSize(); index++) {
+        if (grid->getNeighborsX()[index] == -1)
+            countInvalidX++;
+        if (grid->getNeighborsY()[index] == -1)
+            countInvalidY++;
+        if (grid->getNeighborsZ()[index] == -1)
+            countInvalidZ++;
+    }
+    return { countInvalidX, countInvalidY, countInvalidZ };
+}
+
+std::array<int, 3> testFluidNodeNeighbors(SPtr<Grid> grid)
+{
+    auto countInvalidX = 0;
+    auto countInvalidXY = 0;
+    auto countInvalidXYZ = 0;
+    for (uint index = 0; index < grid->getSize(); index++) {
+        if (grid->getFieldEntry(index) != vf::gpu::FLUID) {
+            continue;
+        }
+
+        auto neighX = grid->getNeighborsX()[index];
+        if (neighX == -1) {
+            countInvalidX++;
+            continue;
+        }
+
+        auto neighXY = grid->getNeighborsY()[neighX];
+        if (neighXY == -1) {
+            countInvalidXY++;
+            continue;
+        }
+
+        auto neighXYZ = grid->getNeighborsZ()[neighXY];
+        if (neighXYZ == -1) {
+            countInvalidXYZ++;
+            continue;
+        }
+    }
+
+    return { countInvalidX, countInvalidXY, countInvalidXYZ };
+}
+
+class findNeighborsIntegrationTest : public ::testing::Test
+{
+protected:
+    SPtr<MultipleGridBuilder> gridBuilder;
+    void SetUp() override
+    {
+        auto gridFactory = GridFactory::make();
+        gridFactory->setTriangularMeshDiscretizationMethod(TriangularMeshDiscretizationMethod::POINT_IN_OBJECT);
+        gridBuilder = MultipleGridBuilder::makeShared(gridFactory);
+
+        // init logger to avoid segmentation fault in buildGrids
+        logging::Logger::addStream(&std::cout);
+        logging::Logger::setDebugLevel(logging::Logger::Level::WARNING);
+        logging::Logger::timeStamp(logging::Logger::ENABLE);
+        logging::Logger::enablePrintedRankNumbers(logging::Logger::ENABLE);
+    }
+};
+
+TEST_F(findNeighborsIntegrationTest, grid1)
+{
+    const real dx = 0.15;
+    gridBuilder->addCoarseGrid(0.0, 0.0, 0.0, 1.0, 1.0, 1.0, dx);
+
+    gridBuilder->buildGrids(LBM, false);
+    auto grid = gridBuilder->getGrid(0);
+
+    // Only the last layer of nodes should have invalid neighbors. The grid is a cube with a side length of 9 nodes
+    // -> 9 * 9 = 81 invalid nodes are expected
+    auto numberOfInvalidNeighbors = countInvalidNeighbors(grid);
+    auto expected = 9 * 9;
+    EXPECT_THAT(numberOfInvalidNeighbors[0], testing::Eq(expected));
+    EXPECT_THAT(numberOfInvalidNeighbors[1], testing::Eq(expected));
+    EXPECT_THAT(numberOfInvalidNeighbors[2], testing::Eq(expected));
+
+    // additional test: all fluid nodes should have valid neighbors
+    auto numberInvalidFluidNeighbors = testFluidNodeNeighbors(grid);
+    EXPECT_THAT(numberInvalidFluidNeighbors[0], testing::Eq(0));
+    EXPECT_THAT(numberInvalidFluidNeighbors[1], testing::Eq(0));
+    EXPECT_THAT(numberInvalidFluidNeighbors[2], testing::Eq(0));
+}
+
+TEST_F(findNeighborsIntegrationTest, grid2)
+{
+    const real dx = 1.0 / 64;
+    gridBuilder->addCoarseGrid(-0.6, -0.6, -0.6, 0.6, 0.6, 0.6, dx);
+
+    gridBuilder->buildGrids(LBM, false);
+    auto grid = gridBuilder->getGrid(0);
+
+    // Only the last layer of nodes should have invalid neighbors. The grid is a cube with a side length of 79 nodes
+    // -> 79 * 79 invalid nodes are expected
+    auto numberOfInvalidNeighbors = countInvalidNeighbors(grid);
+    auto expected = 79 * 79;
+    EXPECT_THAT(numberOfInvalidNeighbors[0], testing::Eq(expected));
+    EXPECT_THAT(numberOfInvalidNeighbors[1], testing::Eq(expected));
+    EXPECT_THAT(numberOfInvalidNeighbors[2], testing::Eq(expected));
+
+    // additional test: all fluid nodes should have valid neighbors
+    auto numberInvalidFluidNeighbors = testFluidNodeNeighbors(grid);
+    EXPECT_THAT(numberInvalidFluidNeighbors[0], testing::Eq(0));
+    EXPECT_THAT(numberInvalidFluidNeighbors[1], testing::Eq(0));
+    EXPECT_THAT(numberInvalidFluidNeighbors[2], testing::Eq(0));
+}
+
+TEST_F(findNeighborsIntegrationTest, validFluidNeighbors1)
+{
+    real dx = 0.17;
+    gridBuilder->addCoarseGrid(0.0, 0.0, 0.0, 1.0, 1.0, 1.0, dx);
+
+    gridBuilder->buildGrids(LBM, false);
+    auto grid = gridBuilder->getGrid(0);
+
+    auto numberInvalidFluidNeighbors = testFluidNodeNeighbors(grid);
+    EXPECT_THAT(numberInvalidFluidNeighbors[0], testing::Eq(0));
+    EXPECT_THAT(numberInvalidFluidNeighbors[1], testing::Eq(0));
+    EXPECT_THAT(numberInvalidFluidNeighbors[2], testing::Eq(0));
+}
+
+TEST_F(findNeighborsIntegrationTest, validFluidNeighbors2)
+{
+    real dx = 0.18;
+    gridBuilder->addCoarseGrid(0.0, 0.0, 0.0, 1.0, 1.0, 1.0, dx);
+
+    gridBuilder->buildGrids(LBM, false);
+    auto grid = gridBuilder->getGrid(0);
+
+    auto numberInvalidFluidNeighbors = testFluidNodeNeighbors(grid);
+    EXPECT_THAT(numberInvalidFluidNeighbors[0], testing::Eq(0));
+    EXPECT_THAT(numberInvalidFluidNeighbors[1], testing::Eq(0));
+    EXPECT_THAT(numberInvalidFluidNeighbors[2], testing::Eq(0));
+}
diff --git a/src/gpu/GridGenerator/grid/NodeValues.h b/src/gpu/GridGenerator/grid/NodeValues.h
index b8312b0673337d11b4bdf0b8052e89d92ce127ef..f1a948cd9a84d1d2454c37c9a23a25639f598e75 100644
--- a/src/gpu/GridGenerator/grid/NodeValues.h
+++ b/src/gpu/GridGenerator/grid/NodeValues.h
@@ -33,9 +33,7 @@
 #ifndef NodeValues_H
 #define NodeValues_H
 
-namespace vf
-{
-namespace gpu
+namespace vf::gpu
 {
 
 static constexpr char FLUID = 0;
@@ -73,7 +71,6 @@ static constexpr char Q_DEPRECATED              = 52;
 
 static constexpr char OVERLAP_TMP = 60;
 
-} // namespace gpu
-} // namespace vf
+}
 
 #endif
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index 92e6883bfe5f4fc2fcab4dbbd98c164a38bf7192..1d2d5677e7014277659fd12e0e2348a5894a7845 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -20,6 +20,7 @@
 #include "Output/AnalysisData.hpp"
 #include "Output/InterfaceDebugWriter.hpp"
 #include "Output/EdgeNodeDebugWriter.hpp"
+#include "Output/NeighborDebugWriter.hpp"
 #include "Output/VeloASCIIWriter.hpp"
 //////////////////////////////////////////////////////////////////////////
 #include "Utilities/Buffer2D.hpp"
@@ -399,6 +400,8 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     // 1000000.0 << " MB\n" << std::endl;
     //////////////////////////////////////////////////////////////////////////
 
+    // NeighborDebugWriter::writeNeighborLinkLinesDebug(para.get());
+
     // InterfaceDebugWriter::writeInterfaceLinesDebugCF(para.get());
     // InterfaceDebugWriter::writeInterfaceLinesDebugFC(para.get());
 
diff --git a/src/gpu/VirtualFluids_GPU/Output/InterfaceDebugWriter.hpp b/src/gpu/VirtualFluids_GPU/Output/InterfaceDebugWriter.hpp
index eb00b43acde44d8e5b3af43843c565b1774e65e6..0b1e9dc1c25457457eabe3013a288c4c93577dc3 100644
--- a/src/gpu/VirtualFluids_GPU/Output/InterfaceDebugWriter.hpp
+++ b/src/gpu/VirtualFluids_GPU/Output/InterfaceDebugWriter.hpp
@@ -3,8 +3,6 @@
 
 #include <fstream>
 #include <sstream>
-#include <stdio.h>
-// #include <math.h>
 #include "Core/StringUtilities/StringUtil.h"
 #include "lbm/constants/D3Q27.h"
 #include "LBM/LB.h"
@@ -906,6 +904,6 @@ void writeRecvNodesStream(Parameter *para)
         WbWriterVtkXmlBinary::getInstance()->writeNodesWithNodeData(filenameVec, nodesVec, datanames, nodedata);
     }
 }
-} // namespace InterfaceDebugWriter
 
+} // namespace InterfaceDebugWriter
 #endif
diff --git a/src/gpu/VirtualFluids_GPU/Output/NeighborDebugWriter.hpp b/src/gpu/VirtualFluids_GPU/Output/NeighborDebugWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..648c1fb5d5480c12f703f7dfbbd717ec1c515da2
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/Output/NeighborDebugWriter.hpp
@@ -0,0 +1,61 @@
+#ifndef NEIGHBORDEBUG_HPP
+#define NEIGHBORDEBUG_HPP
+
+#include "LBM/LB.h"
+#include "Logger.h"
+#include "Parameter/Parameter.h"
+#include "basics/utilities/UbSystem.h"
+#include "grid/NodeValues.h"
+#include "lbm/constants/D3Q27.h"
+#include <basics/writer/WbWriterVtkXmlBinary.h>
+
+#include "Utilities/FindNeighbors.h"
+#include "VirtualFluids_GPU/Communication/Communicator.h"
+#include "Core/StringUtilities/StringUtil.h"
+
+namespace NeighborDebugWriter
+{
+
+inline void writeNeighborLinkLines(Parameter *para, const int level, const uint numberOfNodes, const int direction,
+                                   const std::string &name)
+{
+    VF_LOG_INFO("Write node links in direction {}.", direction);
+    std::vector<UbTupleFloat3> nodes(numberOfNodes * 2);
+    std::vector<UbTupleInt2> cells(numberOfNodes);
+
+    for (uint position = 0; position < numberOfNodes; position++) {
+        if (para->getParH(level)->typeOfGridNode[position] != GEO_FLUID)
+            continue;
+
+        const double x1 = para->getParH(level)->coordinateX[position];
+        const double x2 = para->getParH(level)->coordinateY[position];
+        const double x3 = para->getParH(level)->coordinateZ[position];
+
+        const uint positionNeighbor = getNeighborIndex(para->getParH(level).get(), position, direction);
+
+        const double x1Neighbor = para->getParH(level)->coordinateX[positionNeighbor];
+        const double x2Neighbor = para->getParH(level)->coordinateY[positionNeighbor];
+        const double x3Neighbor = para->getParH(level)->coordinateZ[positionNeighbor];
+
+        nodes.emplace_back(float(x1), float(x2), float(x3));
+        nodes.emplace_back(float(x1Neighbor), float(x2Neighbor), float(x3Neighbor));
+
+        cells.emplace_back((int)nodes.size() - 2, (int)nodes.size() - 1);
+    }
+    WbWriterVtkXmlBinary::getInstance()->writeLines(name, nodes, cells);
+}
+
+inline void writeNeighborLinkLinesDebug(Parameter *para)
+{
+    for (int level = 0; level <= para->getMaxLevel(); level++) {
+        for (int direction = vf::lbm::dir::STARTDIR; direction <= vf::lbm::dir::ENDDIR; direction++) {
+            const std::string fileName = para->getFName() + "_" + StringUtil::toString<int>(level) + "_Link_" +
+                                         std::to_string(direction) + "_Debug.vtk";
+            writeNeighborLinkLines(para, level, para->getParH(level)->numberOfNodes, direction, fileName);
+        }
+    }
+}
+
+} // namespace NeighborDebugWriter
+
+#endif
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
index 87f8533f6c6ce6d40587f55e83aea9a33b09f5ad..709fa092e70d9c3f0da80ab08c8ae68124cacf26 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -663,9 +663,9 @@ void Parameter::setD3Qxx(int d3qxx)
 {
     this->D3Qxx = d3qxx;
 }
-void Parameter::setMaxLevel(int maxlevel)
+void Parameter::setMaxLevel(int numberOfLevels)
 {
-    this->maxlevel = maxlevel - 1;
+    this->maxlevel = numberOfLevels - 1;
     this->fine = this->maxlevel;
     parH.resize(this->maxlevel + 1);
     parD.resize(this->maxlevel + 1);
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
index 435120d378de744ab017c337d822680c71eda43f..3ae501d12c89a9432fe44ad76a1ac329c324bf8a 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
@@ -405,7 +405,7 @@ public:
     void setDiffMod(int DiffMod);
     void setDiffusivity(real Diffusivity);
     void setD3Qxx(int d3qxx);
-    void setMaxLevel(int maxlevel);
+    void setMaxLevel(int numberOfLevels);
     void setParticleBasicLevel(int pbl);
     void setParticleInitLevel(int pil);
     void setNumberOfParticles(int nop);
diff --git a/src/gpu/VirtualFluids_GPU/Utilities/FindNeighbors.h b/src/gpu/VirtualFluids_GPU/Utilities/FindNeighbors.h
new file mode 100644
index 0000000000000000000000000000000000000000..2ecbed5fd6c4b84b92c05953a4c96ecdc3b988de
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/Utilities/FindNeighbors.h
@@ -0,0 +1,33 @@
+#ifndef FIND_NEIGHBORS_H
+#define FIND_NEIGHBORS_H
+
+#include "Parameter/Parameter.h"
+#include "lbm/constants/D3Q27.h"
+
+using namespace vf::lbm::dir;
+
+// Only use for fluid nodes!
+inline uint getNeighborIndex(LBMSimulationParameter *parH, const uint position, const int direction)
+{
+    uint nodeIndex = position;
+
+    if (mapForPointerChasing.at(direction).counterInverse != 0) {
+        nodeIndex = parH->neighborInverse[nodeIndex];
+    }
+
+    for (uint x = 0; x < mapForPointerChasing.at(direction).counterX; x++) {
+        nodeIndex = parH->neighborX[nodeIndex];
+    }
+
+    for (uint y = 0; y < mapForPointerChasing.at(direction).counterY; y++) {
+        nodeIndex = parH->neighborY[nodeIndex];
+    }
+
+    for (uint z = 0; z < mapForPointerChasing.at(direction).counterZ; z++) {
+        nodeIndex = parH->neighborZ[nodeIndex];
+    }
+
+    return nodeIndex;
+}
+
+#endif
diff --git a/src/lbm/constants/D3Q27.h b/src/lbm/constants/D3Q27.h
index 44de7073729d85a802142ea0f9331ec571ff6bb7..4cd9a5be8f070e7eaa18c3114075b4d728f3a73a 100644
--- a/src/lbm/constants/D3Q27.h
+++ b/src/lbm/constants/D3Q27.h
@@ -1,6 +1,9 @@
 #ifndef LBM_D3Q27_H
 #define LBM_D3Q27_H
 
+#include <map>
+#include "basics/Core/DataTypes.h"
+
 namespace vf::lbm::dir
 {
 
@@ -36,7 +39,49 @@ static constexpr int DIR_MMP = 22;	 // TSW
 static constexpr int DIR_PPM = 23;	 // BNE
 static constexpr int DIR_MPM = 24;	 // BNW
 static constexpr int DIR_PMM = 25;	 // BSE
-static constexpr int DIR_MMM = 26;	 // BSW 
+static constexpr int DIR_MMM = 26;	 // BSW
+
+struct countersForPointerChasing{
+    uint counterInverse;
+    uint counterX;
+    uint counterY;
+    uint counterZ;
+};
+
+const std::map<const int, const countersForPointerChasing> mapForPointerChasing = 
+{
+    {DIR_000, countersForPointerChasing{0, 0, 0, 0}},
+    {DIR_P00, countersForPointerChasing{0, 1, 0, 0}},
+    {DIR_M00, countersForPointerChasing{1, 0, 1, 1}},
+    {DIR_0P0, countersForPointerChasing{0, 0, 1, 0}},
+    {DIR_0M0, countersForPointerChasing{1, 1, 0, 1}},
+    {DIR_00P, countersForPointerChasing{0, 0, 0, 1}},
+    {DIR_00M, countersForPointerChasing{1, 1, 1, 0}},
+
+    {DIR_PP0, countersForPointerChasing{0, 1, 1, 0}},
+    {DIR_MM0, countersForPointerChasing{1, 0, 0, 1}},
+    {DIR_PM0, countersForPointerChasing{1, 2, 0, 1}},
+    {DIR_MP0, countersForPointerChasing{1, 0, 2, 1}},
+    {DIR_P0P, countersForPointerChasing{0, 1, 0, 1}},
+    {DIR_M0M, countersForPointerChasing{1, 0, 1, 0}},
+    {DIR_P0M, countersForPointerChasing{1, 2, 1, 0}},
+    {DIR_M0P, countersForPointerChasing{1, 0, 1, 2}},
+    {DIR_0PP, countersForPointerChasing{0, 0, 1, 1}},
+    {DIR_0MM, countersForPointerChasing{1, 1, 0, 0}},
+    {DIR_0PM, countersForPointerChasing{1, 1, 2, 0}},
+    {DIR_0MP, countersForPointerChasing{1, 1, 0, 2}},
+
+    {DIR_PPP, countersForPointerChasing{0, 1, 1, 1}},
+    {DIR_MPP, countersForPointerChasing{1, 0, 2, 2}},
+    {DIR_PMP, countersForPointerChasing{1, 2, 0, 2}},
+    {DIR_MMP, countersForPointerChasing{1, 0, 0, 2}},
+    {DIR_PPM, countersForPointerChasing{1, 2, 2, 0}},
+    {DIR_MPM, countersForPointerChasing{1, 0, 2, 0}},
+    {DIR_PMM, countersForPointerChasing{1, 2, 0, 0}},
+    {DIR_MMM, countersForPointerChasing{1, 0, 0, 0}}
+};
+
+
 
 // used in the CPU version
 // static constexpr int INV_P00 = DIR_M00;