     thrust::device_vector<Real3> sortedDerivTauXyXzYz(numObjectsH->numAllMarkers);
     sortedXSPHandShift.resize(numObjectsH->numAllMarkers);
 
-    thrust::fill(sortedDerivVelRho.begin(), sortedDerivVelRho.end(), mR4(0));
-    thrust::fill(sortedXSPHandShift.begin(), sortedXSPHandShift.end(), mR3(0));
-    if(paramsH->elastic_SPH){
-        thrust::fill(sortedDerivTauXxYyZz.begin(), sortedDerivTauXxYyZz.end(), mR3(0));
-        thrust::fill(sortedDerivTauXyXzYz.begin(), sortedDerivTauXyXzYz.end(), mR3(0));
-    }
-
-    // Find the index which is related to the wall boundary particle
-    uint numBlocks1, numThreads1;
-    computeGridSize((int)numObjectsH->numAllMarkers, 256, numBlocks1, numThreads1);
-    thrust::device_vector<uint> indexOfIndex(numObjectsH->numAllMarkers);
-    thrust::device_vector<uint> identityOfIndex(numObjectsH->numAllMarkers);
-    calIndexOfIndex<<<numBlocks1, numThreads1>>>(U1CAST(indexOfIndex), U1CAST(identityOfIndex), 
-        U1CAST(markersProximityD->gridMarkerIndexD), numObjectsH->numFluidMarkers, 
-        numObjectsH->numBoundaryMarkers, numObjectsH->numAllMarkers);
-    thrust::remove_if(indexOfIndex.begin(), indexOfIndex.end(), identityOfIndex.begin(), 
-        thrust::identity<int>());
-
-    if (density_initialization == 0){
+    // Re-Initialize the density after several time steps
+    if (density_initialization >= paramsH->densityReinit){
         thrust::device_vector<Real4> rhoPresMuD_old = sortedSphMarkersD->rhoPresMuD;
         printf("Re-initializing density after %d steps.\n", paramsH->densityReinit);
-        calcRho_kernel<<<numBlocks, numThreads>>>(U1CAST(indexOfIndex), 
+        calcRho_kernel<<<numBlocks, numThreads>>>(
             mR4CAST(sortedSphMarkersD->posRadD), mR4CAST(sortedSphMarkersD->rhoPresMuD),
             mR4CAST(rhoPresMuD_old), U1CAST(markersProximityD->cellStartD),
             U1CAST(markersProximityD->cellEndD), numObjectsH->numFluidMarkers,
             numObjectsH->numBoundaryMarkers, numObjectsH->numRigid_SphMarkers,
             density_initialization, isErrorD);
-            ChUtilsDevice::Sync_CheckError(isErrorH, isErrorD, "calcRho_kernel");
-    }
+        ChUtilsDevice::Sync_CheckError(isErrorH, isErrorD, "calcRho_kernel");
+        density_initialization = 0;
+    }
+    density_initialization++;
     
+    // Execute the kernel
     if(paramsH->elastic_SPH){ // For granular material
         *isErrorH = false;
         cudaMemcpy(isErrorD, isErrorH, sizeof(bool), cudaMemcpyHostToDevice);
 
         // execute the kernel Navier_Stokes and Shear_Stress_Rate in one kernel
-        NS_SSR<<<numBlocks, numThreads>>>( U1CAST(indexOfIndex), U1CAST(fsiGeneralData->activityIdentifierD),
+        NS_SSR<<<numBlocks, numThreads>>>(U1CAST(fsiGeneralData->activityIdentifierD),
             mR4CAST(sortedDerivVelRho), mR3CAST(sortedDerivTauXxYyZz), mR3CAST(sortedDerivTauXyXzYz),
             mR3CAST(sortedXSPHandShift), mR4CAST(sortedSphMarkersD->posRadD),
             mR3CAST(sortedSphMarkersD->velMasD), mR4CAST(sortedSphMarkersD->rhoPresMuD),
             mR3CAST(bceWorker->velMas_ModifiedBCE), mR4CAST(bceWorker->rhoPreMu_ModifiedBCE),
             mR3CAST(sortedSphMarkersD->tauXxYyZzD), mR3CAST(sortedSphMarkersD->tauXyXzYzD),
             U1CAST(markersProximityD->gridMarkerIndexD), U1CAST(markersProximityD->cellStartD),
-            U1CAST(markersProximityD->cellEndD), numObjectsH->numFluidMarkers,
-            numObjectsH->numBoundaryMarkers, numObjectsH->numRigid_SphMarkers, isErrorD);
+            U1CAST(markersProximityD->cellEndD), U1CAST(markersProximityD->mapOriginalToSorted),
+            numObjectsH->numAllMarkers, isErrorD);
         ChUtilsDevice::Sync_CheckError(isErrorH, isErrorD, "Navier_Stokes and Shear_Stress_Rate");
     }
     else{ // For fluid 
         *isErrorH = false;
         cudaMemcpy(isErrorD, isErrorH, sizeof(bool), cudaMemcpyHostToDevice);
 
+        // Find the index which is related to the wall boundary particle
+        uint numBlocks1, numThreads1;
+        computeGridSize((int)numObjectsH->numFluidMarkers +
+            (int)numObjectsH->numRigid_SphMarkers, 256, numBlocks1, numThreads1);
+
+        thrust::device_vector<uint> indexOfIndex(numObjectsH->numAllMarkers);
+        thrust::device_vector<uint> identityOfIndex(numObjectsH->numAllMarkers);
+        calIndexOfIndex<<<numBlocks, numThreads>>>(U1CAST(indexOfIndex), U1CAST(identityOfIndex), 
+            U1CAST(markersProximityD->gridMarkerIndexD), numObjectsH->numFluidMarkers, 
+            numObjectsH->numBoundaryMarkers, numObjectsH->numAllMarkers);
+        thrust::remove_if(indexOfIndex.begin(), indexOfIndex.end(), identityOfIndex.begin(), 
+            thrust::identity<int>());
+
         // execute the kernel
-        Navier_Stokes<<<numBlocks, numThreads>>>( U1CAST(indexOfIndex), 
+        Navier_Stokes<<<numBlocks1, numThreads1>>>(U1CAST(indexOfIndex), 
             mR4CAST(sortedDerivVelRho), mR3CAST(sortedXSPHandShift), mR4CAST(sortedSphMarkersD->posRadD),
             mR3CAST(sortedSphMarkersD->velMasD), mR4CAST(sortedSphMarkersD->rhoPresMuD),
             mR3CAST(bceWorker->velMas_ModifiedBCE), mR4CAST(bceWorker->rhoPreMu_ModifiedBCE),

}

//--------------------------------------------------------------------------------------------------------------------------------
__global__ void NS_SSR( uint* indexOfIndex,
                        uint* activityIdentifierD,
                        Real4* sortedDerivVelRho,
                        Real3* sortedDerivTauXxYyZz,
                        Real3* sortedDerivTauXyXzYz,

}

//--------------------------------------------------------------------------------------------------------------------------------
__global__ void calcRho_kernel(Real4* sortedPosRad,
                               Real4* sortedRhoPreMu,
                               Real4* sortedRhoPreMu_old,
                               uint* cellStart,
