 }
 
 // -----------------------------------------------------------------------------
-/**
- * @brief ApplyBoundarySPH_Markers
- * @details
- * 		applies periodic boundary conditions in y, and z. The inlet/outlet BC is applied in the x direction.
- * 		This functions needs to be tested.
- */
+// Apply periodic boundary conditions in y, and z. 
+// The inlet/outlet BC is applied in the x direction.
+// This functions needs to be tested.
 void ChFluidDynamics::ApplyModifiedBoundarySPH_Markers(std::shared_ptr<SphMarkerDataD> sphMarkersD) {
     uint numBlocks, numThreads;
     computeGridSize((int)numObjectsH->numAllMarkers, 256, numBlocks, numThreads);
-    ApplyInletBoundaryXKernel<<<numBlocks, numThreads>>>(mR4CAST(sphMarkersD->posRadD), mR3CAST(sphMarkersD->velMasD),
-                                                         mR4CAST(sphMarkersD->rhoPresMuD));
+    ApplyInletBoundaryXKernel<<<numBlocks, numThreads>>>(
+        mR4CAST(sphMarkersD->posRadD), mR3CAST(sphMarkersD->velMasD),
+        mR4CAST(sphMarkersD->rhoPresMuD));
     cudaDeviceSynchronize();
     cudaCheckError();
 
     // these are useful anyway for out of bound particles
-    ApplyPeriodicBoundaryYKernel<<<numBlocks, numThreads>>>(mR4CAST(sphMarkersD->posRadD),
-                                                            mR4CAST(sphMarkersD->rhoPresMuD),
-                                                            U1CAST(fsiSystem.fsiGeneralData->activityIdentifierD));
-    cudaDeviceSynchronize();
-    cudaCheckError();
-
-    ApplyPeriodicBoundaryZKernel<<<numBlocks, numThreads>>>(mR4CAST(sphMarkersD->posRadD),
-                                                            mR4CAST(sphMarkersD->rhoPresMuD),
-                                                            U1CAST(fsiSystem.fsiGeneralData->activityIdentifierD));
+    ApplyPeriodicBoundaryYKernel<<<numBlocks, numThreads>>>(
+        mR4CAST(sphMarkersD->posRadD), mR4CAST(sphMarkersD->rhoPresMuD),
+        U1CAST(fsiSystem.fsiGeneralData->activityIdentifierD));
+    cudaDeviceSynchronize();
+    cudaCheckError();
+
+    ApplyPeriodicBoundaryZKernel<<<numBlocks, numThreads>>>(
+        mR4CAST(sphMarkersD->posRadD), mR4CAST(sphMarkersD->rhoPresMuD),
+        U1CAST(fsiSystem.fsiGeneralData->activityIdentifierD));
     cudaDeviceSynchronize();
     cudaCheckError();
 }

                                  uint* cellStart,
                                  uint* cellEnd,
                                  size_t numAllMarkers) {
    uint index = __mul24(blockIdx.x, blockDim.x) + threadIdx.x;
    if (index >= numAllMarkers)
        return;

    // read particle data from sorted arrays
    Real3 posRadA = mR3(sortedPosRad[index]);
    Real4 rhoPreMuA = sortedRhoPreMu[index];

    // If density initialization should only be applied to fluid particles
    // if (rhoPreMuA.w > -.1)
    //    return;

    // get address in grid
    int3 gridPos = calcGridPos(posRadA);


// -----------------------------------------------------------------------------
// Kernel to apply periodic BC along z
__global__ void ApplyPeriodicBoundaryZKernel(Real4* posRadD, 
                                             Real4* rhoPresMuD, 
                                             uint* activityIdentifierD) {
    uint index = blockIdx.x * blockDim.x + threadIdx.x;
    if (index >= numObjectsD.numAllMarkers)
        return;

    uint activity = activityIdentifierD[index];
    if (activity == 0)
        return; // no need to do anything if it is not an active particle

    Real4 rhoPresMu = rhoPresMuD[index];
    if (fabs(rhoPresMu.w) < .1)
        return; // no need to do anything if it is a boundary particle

    Real3 posRad = mR3(posRadD[index]);
    Real h = posRadD[index].w;

