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

    //------------------------------------------------------------------------
    // thread per particle
    uint numBlocks, numThreads;
    computeGridSize((int)numObjectsH->numFluidMarkers +
        (int)numObjectsH->numRigid_SphMarkers, 256, numBlocks, numThreads);

    // Execute the kernel
    thrust::device_vector<Real4> sortedDerivVelRho(numObjectsH->numAllMarkers);

    thrust::device_vector<Real3> sortedDerivTauXyXzYz(numObjectsH->numAllMarkers);
    sortedXSPHandShift.resize(numObjectsH->numAllMarkers);

    // 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>>>(
            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");
        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(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), 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<<<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),
