diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index 6a9fe44f756c1b4fde106e60953a540a161c8e81..aa9df75c5ec56b6114e7f5d5983a491e6bc39388 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -101,7 +101,7 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     gridProvider.setBoundingBox();
 
     para->setRe(para->getVelocity() * (real)1.0 / para->getViscosity());
-    para->setlimitOfNodesForVTK(30000000); // max 30 Million nodes per VTK file
+    para->setLimitOfNodesForVTK(30000000); // max 30 Million nodes per VTK file
     if (para->getDoRestart())
         para->setStartTurn(para->getTimeDoRestart());
     else
diff --git a/src/gpu/VirtualFluids_GPU/Output/FileWriter.cpp b/src/gpu/VirtualFluids_GPU/Output/FileWriter.cpp
index bb0f1b83edf27be2964dcc46f91fff2e61b32659..97f712aeb97bc901b271eeddaafda499457fca3e 100644
--- a/src/gpu/VirtualFluids_GPU/Output/FileWriter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Output/FileWriter.cpp
@@ -284,7 +284,7 @@ std::vector<std::string> FileWriter::writeUnstructuredGridLT(std::shared_ptr<Par
 
                 WriterUtilities::getIndicesOfAllNodesInOct(indicesOfOct, pos, para->getParHConst(level).get());
 
-                if (WriterUtilities::isPeriodicCell(para.get(), level, indicesOfOct[0], indicesOfOct[6])){
+                if (WriterUtilities::isPeriodicCell(para.get(), level, indicesOfOct[0], indicesOfOct[6])) {
                     continue;
                 }
 
@@ -293,7 +293,7 @@ std::vector<std::string> FileWriter::writeUnstructuredGridLT(std::shared_ptr<Par
 
                 //////////////////////////////////////////////////////////////////////////
                 if (allNodesValid) {
-                    WriterUtilities::calculateRelativePositionInPart(relativePosInPart, indicesOfOct, startPosition);
+                    WriterUtilities::calculateRelativeNodeIndexInPart(relativePosInPart, indicesOfOct, startPosition);
                     cells.push_back(makeUbTupleFromArray(relativePosInPart));
                 }
             }
diff --git a/src/gpu/VirtualFluids_GPU/Output/WriterUtilities.cpp b/src/gpu/VirtualFluids_GPU/Output/WriterUtilities.cpp
index 50ad3c08ad3289631ccc5f2f279cf9aaa36e1190..d6775589866843d45c306fbe00ccf20561b45967 100644
--- a/src/gpu/VirtualFluids_GPU/Output/WriterUtilities.cpp
+++ b/src/gpu/VirtualFluids_GPU/Output/WriterUtilities.cpp
@@ -13,40 +13,43 @@ uint WriterUtilities::calculateNumberOfParts(const Parameter* parameter, uint le
     return (uint)parameter->getParHConst(level)->numberOfNodes / parameter->getLimitOfNodesForVTK() + 1;
 }
 
-bool WriterUtilities::isPeriodicCell(const Parameter* para, int level, unsigned int number1, unsigned int number7)
+bool WriterUtilities::isPeriodicCell(const Parameter* para, int level, unsigned int baseNodeOfCell,
+                                     unsigned int otherNodeInCell)
 {
-    real distance =
-        sqrt(pow(para->getParHConst(level)->coordinateX[number7] - para->getParHConst(level)->coordinateX[number1], 2.) +
-             pow(para->getParHConst(level)->coordinateY[number7] - para->getParHConst(level)->coordinateY[number1], 2.) +
-             pow(para->getParHConst(level)->coordinateZ[number7] - para->getParHConst(level)->coordinateZ[number1], 2.));
-    return distance > 1.01 * sqrt(3 * para->getParHConst(level)->gridSpacing);
+    real distance = sqrt(
+        pow(para->getParHConst(level)->coordinateX[otherNodeInCell] - para->getParHConst(level)->coordinateX[baseNodeOfCell], 2.) +
+        pow(para->getParHConst(level)->coordinateY[otherNodeInCell] - para->getParHConst(level)->coordinateY[baseNodeOfCell], 2.) +
+        pow(para->getParHConst(level)->coordinateZ[otherNodeInCell] - para->getParHConst(level)->coordinateZ[baseNodeOfCell], 2.));
+    return distance > 1.01 * sqrt(3 * pow(para->getParHConst(level)->gridSpacing, 2.));
 }
 
 uint WriterUtilities::calculateNumberOfNodesInPart(const Parameter* para, uint level, uint part)
 {
+    if (part >= WriterUtilities::calculateNumberOfParts(para, level))
+        throw std::runtime_error("The number of nodes for a non-existing part can not be calculated");
     if (((part + 1) * para->getLimitOfNodesForVTK()) > (uint)para->getParHConst(level)->numberOfNodes)
         return (uint)para->getParHConst(level)->numberOfNodes - (part * para->getLimitOfNodesForVTK());
     return para->getLimitOfNodesForVTK();
 }
 
-void WriterUtilities::getIndicesOfAllNodesInOct(std::array<uint, 8>& indices, uint baseNodeOfOct,
+void WriterUtilities::getIndicesOfAllNodesInOct(std::array<uint, 8>& nodeIndices, uint baseNodeOfOct,
                                                 const LBMSimulationParameter* parH)
 {
-    indices[0] = baseNodeOfOct;
-    indices[1] = parH->neighborX[indices[0]];
-    indices[2] = parH->neighborY[indices[1]];
-    indices[3] = parH->neighborY[indices[0]];
-    indices[4] = parH->neighborZ[indices[0]];
-    indices[5] = parH->neighborZ[indices[1]];
-    indices[6] = parH->neighborZ[indices[2]];
-    indices[7] = parH->neighborZ[indices[3]];
+    nodeIndices[0] = baseNodeOfOct;
+    nodeIndices[1] = parH->neighborX[nodeIndices[0]];
+    nodeIndices[2] = parH->neighborY[nodeIndices[1]];
+    nodeIndices[3] = parH->neighborY[nodeIndices[0]];
+    nodeIndices[4] = parH->neighborZ[nodeIndices[0]];
+    nodeIndices[5] = parH->neighborZ[nodeIndices[1]];
+    nodeIndices[6] = parH->neighborZ[nodeIndices[2]];
+    nodeIndices[7] = parH->neighborZ[nodeIndices[3]];
 }
 
-void WriterUtilities::calculateRelativePositionInPart(std::array<uint, 8>& relativePositionInPart,
-                                                      const std::array<uint, 8>& indicesOfOct, uint startPosition)
+void WriterUtilities::calculateRelativeNodeIndexInPart(std::array<uint, 8>& relativePositionInPart,
+                                                      const std::array<uint, 8>& indicesOfOct, uint startPositionOfPart)
 {
     for (size_t i = 0; i < relativePositionInPart.size(); i++) {
-        relativePositionInPart[i] = indicesOfOct[i] - startPosition;
+        relativePositionInPart[i] = indicesOfOct[i] - startPositionOfPart;
     }
 }
 
diff --git a/src/gpu/VirtualFluids_GPU/Output/WriterUtilities.h b/src/gpu/VirtualFluids_GPU/Output/WriterUtilities.h
index 90081fa975c227c6c1538f3d441f80a7c64ed561..35705e4c9414c7c0f4fb43e383364de1329b8d0f 100644
--- a/src/gpu/VirtualFluids_GPU/Output/WriterUtilities.h
+++ b/src/gpu/VirtualFluids_GPU/Output/WriterUtilities.h
@@ -13,14 +13,15 @@ class WriterUtilities
 public:
     static uint calculateNumberOfParts(const Parameter* para, uint level);
     static uint calculateNumberOfNodesInPart(const Parameter* para, uint level, uint part);
-    static std::string makePartFileNameEnding(uint level, int processID, int part, int timestep);
-    static void getIndicesOfAllNodesInOct(std::array<uint, 8>& indices, uint baseNodeOfOct,
+    static bool isPeriodicCell(const Parameter* para, int level, unsigned int baseNodeOfCell, unsigned int otherNodeInCell);
+    static void getIndicesOfAllNodesInOct(std::array<uint, 8>& nodeIndices, uint baseNodeOfOct,
                                           const LBMSimulationParameter* parH);
-    static void calculateRelativePositionInPart(std::array<uint, 8>& relativePositionInPart,
-                                                const std::array<uint, 8>& indicesOfOct, uint startPosition);
+    static void calculateRelativeNodeIndexInPart(std::array<uint, 8>& relativePositionInPart,
+                                                 const std::array<uint, 8>& indicesOfOct, uint startPositionOfPart);
     static bool areAllNodesInOctValidForWriting(const std::array<uint, 8>& indicesOfOct, const LBMSimulationParameter* parH,
                                                 uint endPositionOfPart);
-    static bool isPeriodicCell(const Parameter* para, int level, unsigned int number1, unsigned int number7);
+
+    static std::string makePartFileNameEnding(uint level, int processID, int part, int timestep);
 };
 
 #endif
diff --git a/src/gpu/VirtualFluids_GPU/Output/WriterUtilitiesTest.cpp b/src/gpu/VirtualFluids_GPU/Output/WriterUtilitiesTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7702c5f33fdb315e2ee332629cb775b800297302
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/Output/WriterUtilitiesTest.cpp
@@ -0,0 +1,252 @@
+
+#include "WriterUtilities.h"
+#include "gpu/VirtualFluids_GPU/Utilities/testUtilitiesGPU.h"
+#include <gmock/gmock.h>
+#include <numeric>
+
+TEST(WriterUtilitiesTest, calculateNumberOfParts)
+{
+    const uint level = 0;
+    auto parameter = testingVF::createParameterForLevel(level);
+    parameter->setLimitOfNodesForVTK(10);
+
+    parameter->getParH(level)->numberOfNodes = 23;
+    EXPECT_THAT(WriterUtilities::calculateNumberOfParts(parameter.get(), level), testing::Eq(3));
+
+    parameter->getParH(level)->numberOfNodes = 13;
+    EXPECT_THAT(WriterUtilities::calculateNumberOfParts(parameter.get(), level), testing::Eq(2));
+
+    parameter->getParH(level)->numberOfNodes = 3;
+    EXPECT_THAT(WriterUtilities::calculateNumberOfParts(parameter.get(), level), testing::Eq(1));
+}
+
+TEST(WriterUtilitiesTest, calculateNumberOfNodesInPart)
+{
+    const uint level = 0;
+    auto parameter = testingVF::createParameterForLevel(level);
+    parameter->getParH(level)->numberOfNodes = 13;
+    parameter->setLimitOfNodesForVTK(10);
+
+    uint part = 0;
+    EXPECT_THAT(WriterUtilities::calculateNumberOfNodesInPart(parameter.get(), level, part), testing::Eq(10));
+
+    part = 1;
+    EXPECT_THAT(WriterUtilities::calculateNumberOfNodesInPart(parameter.get(), level, part), testing::Eq(3));
+
+    part = 2;
+    EXPECT_THROW(WriterUtilities::calculateNumberOfNodesInPart(parameter.get(), level, part), std::runtime_error);
+}
+
+class WriterUtilitiesPeriodicCellTest : public testing::Test
+{
+protected:
+    SPtr<Parameter> parameter;
+    SPtr<LBMSimulationParameter> parH;
+    const uint level = 0;
+    const uint baseNodeIndex = 0;
+    const uint otherNodeIndex = 1;
+    std::array<real, 2> coordinates = {0.0, 1.0};
+
+    void SetUp() override
+    {
+        // create a domain with only three layers of nodes
+        // nodes are at the coordinates 0.0, 1.0 and 2.0
+        parameter = testingVF::createParameterForLevel(level);
+        parH = parameter->getParH(level);
+        parH->gridSpacing = 1.0;
+
+        parH->coordinateX = new real[2];
+        parH->coordinateY = new real[2];
+        parH->coordinateZ = new real[2];
+
+        parH->coordinateX[baseNodeIndex] = coordinates[baseNodeIndex];
+        parH->coordinateY[baseNodeIndex] = coordinates[baseNodeIndex];
+        parH->coordinateZ[baseNodeIndex] = coordinates[baseNodeIndex];
+        parH->coordinateX[otherNodeIndex] = coordinates[otherNodeIndex];
+        parH->coordinateY[otherNodeIndex] = coordinates[otherNodeIndex];
+        parH->coordinateZ[otherNodeIndex] = coordinates[otherNodeIndex];
+    }
+
+    void TearDown() override
+    {
+        delete[] parH->coordinateX;
+        delete[] parH->coordinateY;
+        delete[] parH->coordinateZ;
+    }
+};
+
+TEST_F(WriterUtilitiesPeriodicCellTest, cellIsNotPeriodic)
+{
+    EXPECT_FALSE(WriterUtilities::isPeriodicCell(parameter.get(), level, baseNodeIndex, otherNodeIndex));
+    EXPECT_FALSE(WriterUtilities::isPeriodicCell(parameter.get(), level, otherNodeIndex, baseNodeIndex));
+}
+
+TEST_F(WriterUtilitiesPeriodicCellTest, cellIsPeriodicInX)
+{
+    parH->coordinateX[1] = 2.0;
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, baseNodeIndex, otherNodeIndex));
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, otherNodeIndex, baseNodeIndex));
+}
+
+TEST_F(WriterUtilitiesPeriodicCellTest, cellIsPeriodicInY)
+{
+    parH->coordinateY[1] = 2.0;
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, baseNodeIndex, otherNodeIndex));
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, otherNodeIndex, baseNodeIndex));
+}
+
+TEST_F(WriterUtilitiesPeriodicCellTest, cellIsPeriodicInZ)
+{
+    parH->coordinateZ[1] = 2.0;
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, baseNodeIndex, otherNodeIndex));
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, otherNodeIndex, baseNodeIndex));
+}
+TEST_F(WriterUtilitiesPeriodicCellTest, cellIsPeriodicInXY)
+{
+    parH->coordinateX[1] = 2.0;
+    parH->coordinateY[1] = 2.0;
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, baseNodeIndex, otherNodeIndex));
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, otherNodeIndex, baseNodeIndex));
+}
+
+TEST_F(WriterUtilitiesPeriodicCellTest, cellIsPeriodicInXZ)
+{
+    parH->coordinateX[1] = 2.0;
+    parH->coordinateZ[1] = 2.0;
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, baseNodeIndex, otherNodeIndex));
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, otherNodeIndex, baseNodeIndex));
+}
+
+TEST_F(WriterUtilitiesPeriodicCellTest, cellIsPeriodicInYZ)
+{
+    parH->coordinateY[1] = 2.0;
+    parH->coordinateZ[1] = 2.0;
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, baseNodeIndex, otherNodeIndex));
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, otherNodeIndex, baseNodeIndex));
+}
+
+TEST_F(WriterUtilitiesPeriodicCellTest, cellIsPeriodicInXYZ)
+{
+    parH->coordinateX[1] = 2.0;
+    parH->coordinateY[1] = 2.0;
+    parH->coordinateZ[1] = 2.0;
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, baseNodeIndex, otherNodeIndex));
+    EXPECT_TRUE(WriterUtilities::isPeriodicCell(parameter.get(), level, otherNodeIndex, baseNodeIndex));
+}
+
+class WriterUtilitiesNeighborOctTest : public testing::Test
+{
+
+    static void setUpNeighborsNeighborsForOct(LBMSimulationParameter* parH, const std::array<uint, 8>& nodeIndices)
+    {
+        // node indices: MMM, PMM, PPM, MPM,
+        //               MMP, PMP, PPP, MPP
+
+        for (uint i = 0; i < (uint)nodeIndices.size(); i++) {
+            const uint currentNodeIndex = nodeIndices[i];
+            if (i < 4)
+                parH->neighborZ[currentNodeIndex] = nodeIndices[i + 4];
+            else
+                parH->neighborZ[currentNodeIndex] = 99;
+
+            if (i == 0 || i == 4)
+                parH->neighborY[currentNodeIndex] = nodeIndices[i + 3];
+            else if (i == 1 || i == 5)
+                parH->neighborY[currentNodeIndex] = nodeIndices[i + 1];
+            else
+                parH->neighborY[currentNodeIndex] = 999;
+
+            if (i == 0 || i == 4)
+                parH->neighborX[currentNodeIndex] = nodeIndices[i + 1];
+            else if (i == 3 || i == 7)
+                parH->neighborX[currentNodeIndex] = nodeIndices[i - 1];
+            else
+                parH->neighborX[currentNodeIndex] = 9999;
+        }
+    }
+
+public:
+    std::unique_ptr<LBMSimulationParameter> parH = std::make_unique<LBMSimulationParameter>();
+    std::array<uint, 8> nodeIndices;
+
+    void SetUp() override
+    {
+        // set up some node indices from 0 to 7
+        std::iota(nodeIndices.begin(), nodeIndices.end(), 0);
+        std::reverse(nodeIndices.begin(), nodeIndices.end());
+
+        parH->neighborX = new uint[8];
+        parH->neighborY = new uint[8];
+        parH->neighborZ = new uint[8];
+        setUpNeighborsNeighborsForOct(parH.get(), nodeIndices);
+    }
+
+    void TearDown() override
+    {
+        delete[] parH->neighborX;
+        delete[] parH->neighborY;
+        delete[] parH->neighborZ;
+    }
+};
+
+TEST_F(WriterUtilitiesNeighborOctTest, getIndicesOfAllNodesInOct)
+{
+    std::array<uint, 8> resultingNodeIndices;
+    WriterUtilities::getIndicesOfAllNodesInOct(resultingNodeIndices, nodeIndices[0], parH.get());
+    for (uint i = 0; i < 8; i++)
+        EXPECT_THAT(resultingNodeIndices[i], testing::Eq(nodeIndices[i])) << "for index i = " << i << " in nodeIndices";
+}
+
+TEST(WriterUtilitiesTest, calculateRelativeNodeIndexInPart)
+{
+    std::array<uint, 8> indicesOfOct = { 10, 12, 14, 16, 20, 22, 24, 26 };
+    uint startPositionOfPart = 10;
+    std::array<uint, 8> expected = { 0, 2, 4, 6, 10, 12, 14, 16 };
+
+    std::array<uint, 8> result;
+    WriterUtilities::calculateRelativeNodeIndexInPart(result, indicesOfOct, startPositionOfPart);
+    EXPECT_THAT(result, testing::Eq(expected));
+}
+
+class WriterUtilitiesTestNodeValidity : public testing::Test
+{
+protected:
+    std::unique_ptr<LBMSimulationParameter> parH = std::make_unique<LBMSimulationParameter>();
+    std::array<uint, 8> nodeIndices;
+    std::array<uint, 8> typeOfGridNode;
+
+    void SetUp() override
+    {
+        // set up node indices from 0 to 7
+        std::iota(nodeIndices.begin(), nodeIndices.end(), 0);
+
+        std::fill(typeOfGridNode.begin(), typeOfGridNode.end(), GEO_FLUID);
+        parH->typeOfGridNode = typeOfGridNode.data();
+    }
+};
+
+TEST_F(WriterUtilitiesTestNodeValidity, allNodesInOctValidForWriting)
+{
+    uint endPositionOfPart = 7;
+    EXPECT_TRUE(WriterUtilities::areAllNodesInOctValidForWriting(nodeIndices, parH.get(), endPositionOfPart));
+}
+
+TEST_F(WriterUtilitiesTestNodeValidity, areAllNodesInOctValidForWriting_NodeOutOfPart)
+{
+    uint endPositionOfPart = 6;
+    EXPECT_FALSE(WriterUtilities::areAllNodesInOctValidForWriting(nodeIndices, parH.get(), endPositionOfPart));
+}
+
+TEST_F(WriterUtilitiesTestNodeValidity, areAllNodesInOctValidForWriting_NonFluidNode)
+{
+    uint endPositionOfPart = 7;
+    typeOfGridNode[0] = GEO_SOLID;
+    EXPECT_FALSE(WriterUtilities::areAllNodesInOctValidForWriting(nodeIndices, parH.get(), endPositionOfPart));
+}
+
+TEST_F(WriterUtilitiesTestNodeValidity, areAllNodesInOctValidForWriting_NonFluidNodeAtEnd)
+{
+    uint endPositionOfPart = 7;
+    typeOfGridNode[7] = GEO_SOLID;
+    EXPECT_FALSE(WriterUtilities::areAllNodesInOctValidForWriting(nodeIndices, parH.get(), endPositionOfPart));
+}
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
index 4f6a6696b6b6a8660ac046bf7f2af2789606ac86..a3e975d4f4fb93c70fa5978281a9dfdd4fbb2494 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -665,7 +665,7 @@ void Parameter::setOutputCount(unsigned int outputCount)
 {
     this->outputCount = outputCount;
 }
-void Parameter::setlimitOfNodesForVTK(unsigned int limitOfNodesForVTK)
+void Parameter::setLimitOfNodesForVTK(unsigned int limitOfNodesForVTK)
 {
     this->limitOfNodesForVTK = limitOfNodesForVTK;
 }
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
index 2fc77326a9c9c17f8f977ee37bd35e55658a1e08..5193f15e278e870410ca1ac1c9a34fb3904d6743 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
@@ -488,7 +488,7 @@ public:
     void setAngularVelocity(real inAngVel);
     void setStepEnsight(unsigned int step);
     void setOutputCount(unsigned int outputCount);
-    void setlimitOfNodesForVTK(unsigned int limitOfNodesForVTK);
+    void setLimitOfNodesForVTK(unsigned int limitOfNodesForVTK);
     void setStartTurn(unsigned int inStartTurn);
     void setDiffOn(bool isDiff);
     void setCompOn(bool isComp);
diff --git a/src/lbm/constants/D3Q27.h b/src/lbm/constants/D3Q27.h
index 17b663169997bc384ab879a6a6fde6421ec0aee5..8b0c8febd0143540a70541c5c9802f9e606f4b9b 100644
--- a/src/lbm/constants/D3Q27.h
+++ b/src/lbm/constants/D3Q27.h
@@ -10,6 +10,7 @@ namespace vf::lbm::dir
 
 static constexpr size_t STARTDIR = 0;
 static constexpr size_t ENDDIR = 26;
+static constexpr size_t NUMBER_Of_DIRECTIONS = ENDDIR + 1;
 
 // used in the CPU and the GPU version
 static constexpr size_t DIR_000 = 0;