             U1CAST(indexOfIndex), mR4CAST(sortedDerivVelRho), 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),
+            mR4CAST(bceWorker->rhoPreMu_ModifiedBCE), U1CAST(markersProximityD->gridMarkerIndexD),
             U1CAST(markersProximityD->cellStartD), U1CAST(markersProximityD->cellEndD), 
-            numObjectsH->numFluidMarkers, numObjectsH->numBoundaryMarkers, numObjectsH->numRigidMarkers,
-            numObjectsH->numAllMarkers, isErrorD);
+            numObjectsH->numFluidMarkers, numObjectsH->numBoundaryMarkers, 
+            numObjectsH->numRigidMarkers, numObjectsH->numAllMarkers, isErrorD);
         ChUtilsDevice::Sync_CheckError(isErrorH, isErrorD, "Navier_Stokes");
     }
 


    // Viscosity
    Real rAB_Dot_GradWh = dot(dist3, gradW);
    Real rAB_Dot_GradWh_OverDist = rAB_Dot_GradWh / (d * d + paramsD.epsMinMarkersDis * paramsD.HSML * paramsD.HSML);
    Real3 derivV = -paramsD.markerMass *
                       (rhoPresMuA.y / (rhoPresMuA.x * rhoPresMuA.x) + rhoPresMuB.y / (rhoPresMuB.x * rhoPresMuB.x)) *
                       gradW +
                   paramsD.markerMass * (8.0f * multViscosity) * paramsD.mu0 * rAB_Dot_GradWh_OverDist *
                       (velMasA - velMasB) / square(rhoPresMuA.x + rhoPresMuB.x);

    // Artificial viscosity
    Real vAB_Dot_rAB = dot(velMasA - velMasB, dist3);
    if ((vAB_Dot_rAB < 0.0) && (1 == 0)) {  // change to 1==1 if needs artificial viscosity
        Real alpha = paramsD.Ar_vis_alpha;
        Real c_ab = paramsD.Cs;
        Real rho = 0.5f * (rhoPresMuA.x * rhoPresMuB.x);
        Real nu = -alpha * paramsD.HSML * c_ab / rho;
        Real derivM1 =
            -paramsD.markerMass * (nu * vAB_Dot_rAB / (d * d + paramsD.epsMinMarkersDis * paramsD.HSML * paramsD.HSML));
        derivV.x += derivM1 * gradW.x;
        derivV.y += derivM1 * gradW.y;
        derivV.z += derivM1 * gradW.z;

        return;

    uint index = indexOfIndex[id];
    
    // Do nothing for fixed wall BCE particles 
    if (sortedRhoPreMu[index].w > -0.5 && sortedRhoPreMu[index].w < 0.5) {
        sortedDerivVelRho[index] = mR4(0.0);
        return;
