     *isErrorH = false;
     cudaMemcpy(isErrorD, isErrorH, sizeof(bool), cudaMemcpyHostToDevice);
     Update_Fluid_State<<<numBlocks, numThreads>>>(
-        mR3CAST(fsiSystem.fsiGeneralData->vel_XSPH_D), mR3CAST(fsiSystem.fsiGeneralData->vis_vel_SPH_D),
-        mR4CAST(sphMarkersD->posRadD), mR3CAST(sphMarkersD->velMasD), mR4CAST(sphMarkersD->rhoPresMuD), updatePortion,
+        mR3CAST(fsiSystem.fsiGeneralData->vel_XSPH_D), 
+        mR4CAST(sphMarkersD->posRadD), mR3CAST(sphMarkersD->velMasD), 
+        mR4CAST(sphMarkersD->rhoPresMuD), updatePortion,
         numObjectsH->numAllMarkers, paramsH->dT, isErrorD);
     cudaDeviceSynchronize();
     cudaCheckError();
 
     cudaMemcpy(isErrorH, isErrorD, sizeof(bool), cudaMemcpyDeviceToHost);
-    if (*isErrorH == true) {
+    if (*isErrorH == true)
         throw std::runtime_error("Error! program crashed in Update_Fluid_State!\n");
-    }
-    //------------------------------------------------------------------------
-
     cudaFree(isErrorD);
     free(isErrorH);
 }
 
 // -----------------------------------------------------------------------------
-/**
- * @brief ApplyBoundarySPH_Markers
- * @details
- * 		applies periodic boundary conditions in x, y, and z directions
- */
+// Apply periodic boundary conditions in x, y, and z directions
 void ChFluidDynamics::ApplyBoundarySPH_Markers(std::shared_ptr<SphMarkerDataD> sphMarkersD) {
     uint numBlocks, numThreads;
 
     computeGridSize((int)numObjectsH->numAllMarkers, 256, numBlocks, numThreads);
-    ApplyPeriodicBoundaryXKernel<<<numBlocks, numThreads>>>(mR4CAST(sphMarkersD->posRadD),
-                                                            mR4CAST(sphMarkersD->rhoPresMuD),
-                                                            U1CAST(fsiSystem.fsiGeneralData->activityIdentifierD));
-    cudaDeviceSynchronize();
-    cudaCheckError();
-
-    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));
+    ApplyPeriodicBoundaryXKernel<<<numBlocks, numThreads>>>(
+        mR4CAST(sphMarkersD->posRadD), mR4CAST(sphMarkersD->rhoPresMuD),
+        U1CAST(fsiSystem.fsiGeneralData->activityIdentifierD));
+    cudaDeviceSynchronize();
+    cudaCheckError();
+
+    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();
 

    for (uint num = 0; num < numRigidBodies; num++) {
        Real3 detPos = posRadA - posRigidBodiesD[num];
        Real3 Acdomain = paramsD.bodyActiveDomain;
        Real3 ExAcdomain = paramsD.bodyActiveDomain + mR3(2 * RESOLUTION_LENGTH_MULT * paramsD.HSML);
        if (abs(detPos.x) > Acdomain.x || abs(detPos.y) > Acdomain.y || abs(detPos.z) > Acdomain.z)
            isNotActive = isNotActive + 1;
        if (abs(detPos.x) > ExAcdomain.x || abs(detPos.y) > ExAcdomain.y || abs(detPos.z) > ExAcdomain.z)
            isNotExtended = isNotExtended + 1;
    }



// -----------------------------------------------------------------------------
// Kernel to keep particle inside the simulation domain
__global__ void ApplyOutOfBoundaryKernel(Real4* posRadD, 
                                         Real4* rhoPresMuD, 
                                         Real3* velMasD) {
    uint index = blockIdx.x * blockDim.x + threadIdx.x;
    if (index >= numObjectsD.numAllMarkers)
        return;

    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]);
    Real3 vel = mR3(velMasD[index]);
    Real h = posRadD[index].w;
    
    if (posRad.x > 0.5 * paramsD.boxDimX)
        posRad.x = 0.5 * paramsD.boxDimX;
    if (posRad.x < -0.5 * paramsD.boxDimX)
        posRad.x = -0.5 * paramsD.boxDimX;
    if (posRad.y > 0.5 * paramsD.boxDimY)
        posRad.y = 0.5 * paramsD.boxDimY;
    if (posRad.y < -0.5 * paramsD.boxDimY)
        posRad.y = -0.5 * paramsD.boxDimY;
    if (posRad.z > 1.0 * paramsD.boxDimZ)
        posRad.z = 1.0 * paramsD.boxDimZ;
    if (posRad.z < -0.0 * paramsD.boxDimZ)
        posRad.z = -0.0 * paramsD.boxDimZ;

    posRadD[index] = mR4(posRad, h);
    velMasD[index] = mR3(vel);
    return;
