     Real MassOverRho = Mass * paramsD.invrho0 * paramsD.invrho0;
     Real3 MA_gradW = gradW * MassOverRho;
 
-    Real derivVx = (tauXxYyZz_A.x + tauXxYyZz_B.x) * MA_gradW.x + (tauXyXzYz_A.x + tauXyXzYz_B.x) * MA_gradW.y +
-                   (tauXyXzYz_A.y + tauXyXzYz_B.y) * MA_gradW.z;
-    Real derivVy = (tauXyXzYz_A.x + tauXyXzYz_B.x) * MA_gradW.x + (tauXxYyZz_A.y + tauXxYyZz_B.y) * MA_gradW.y +
-                   (tauXyXzYz_A.z + tauXyXzYz_B.z) * MA_gradW.z;
-    Real derivVz = (tauXyXzYz_A.y + tauXyXzYz_B.y) * MA_gradW.x + (tauXyXzYz_A.z + tauXyXzYz_B.z) * MA_gradW.y +
-                   (tauXxYyZz_A.z + tauXxYyZz_B.z) * MA_gradW.z;
+    Real derivVx = (tauXxYyZz_A.x + tauXxYyZz_B.x) * MA_gradW.x + 
+        (tauXyXzYz_A.x + tauXyXzYz_B.x) * MA_gradW.y +
+        (tauXyXzYz_A.y + tauXyXzYz_B.y) * MA_gradW.z;
+    Real derivVy = (tauXyXzYz_A.x + tauXyXzYz_B.x) * MA_gradW.x + 
+        (tauXxYyZz_A.y + tauXxYyZz_B.y) * MA_gradW.y +
+        (tauXyXzYz_A.z + tauXyXzYz_B.z) * MA_gradW.z;
+    Real derivVz = (tauXyXzYz_A.y + tauXyXzYz_B.y) * MA_gradW.x + 
+        (tauXyXzYz_A.z + tauXyXzYz_B.z) * MA_gradW.y +
+        (tauXxYyZz_A.z + tauXxYyZz_B.z) * MA_gradW.z;
 
     // TODO: Visco-plastic model
     // Real vel = length(velMasA);

ChFsiForceExplicitSPH::~ChFsiForceExplicitSPH() {}

//--------------------------------------------------------------------------------------------------------------------------------

void ChFsiForceExplicitSPH::Initialize() {
    ChFsiForce::Initialize();
    cudaMemcpyToSymbolAsync(paramsD, paramsH.get(), sizeof(SimParams));

                        Real3 velMasB = sortedVelMas[j];
                        if (rhoPresMuB.w > -0.5) {
                            int bceIndexB = gridMarkerIndex[j] - numObjectsD.numFluidMarkers;
                            if (!(bceIndexB >= 0 && bceIndexB < 
                                numObjectsD.numAllMarkers - numObjectsD.numFluidMarkers)) {
                                printf("Error! bceIndex out of bound, collideCell !\n");
                            }
                            rhoPresMuB = rhoPreMu_ModifiedBCE[bceIndexB];
                            velMasB = velMas_ModifiedBCE[bceIndexB];
                        }
                        Real multViscosit = 1;

                        derivVelRho += DifVelocityRho(Gi, dist3, d, sortedPosRad[index], sortedPosRad[j],
                             velMasA, velMasB, rhoPresMuA, rhoPresMuB, multViscosit);
                        preGra += GradientOperator(Gi, dist3, sortedPosRad[index], sortedPosRad[j], 
                            -rhoPresMuA.y, rhoPresMuB.y, rhoPresMuA, rhoPresMuB);
                        velxGra += GradientOperator(Gi, dist3, sortedPosRad[index], sortedPosRad[j], 
                            velMasA.x, velMasB.x, rhoPresMuA, rhoPresMuB);
                        velyGra += GradientOperator(Gi, dist3, sortedPosRad[index], sortedPosRad[j], 
                            velMasA.y, velMasB.y, rhoPresMuA, rhoPresMuB);
                        velzGra += GradientOperator(Gi, dist3, sortedPosRad[index], sortedPosRad[j], 
                            velMasA.z, velMasB.z, rhoPresMuA, rhoPresMuB);
                        velxLap += LaplacianOperator(Gi, Li, dist3, sortedPosRad[index], sortedPosRad[j], 
                            velMasA.x, velMasB.x, rhoPresMuA, rhoPresMuB);
                        velyLap += LaplacianOperator(Gi, Li, dist3, sortedPosRad[index], sortedPosRad[j], 
                            velMasA.y, velMasB.y, rhoPresMuA, rhoPresMuB);
                        velzLap += LaplacianOperator(Gi, Li, dist3, sortedPosRad[index], sortedPosRad[j], 
                            velMasA.z, velMasB.z, rhoPresMuA, rhoPresMuB);

                        if (d > paramsD.HSML * 1.0e-9 && sum_w_i < paramsD.C_Wi) {
                            sum_w_i = sum_w_i + W3h(d, sortedPosRad[index].w) * paramsD.volume0;
