 				real3 posRadB = FETCH(sortedPosRad, j);
 				real3 dist3Alpha = posRadA - posRadB;
 				real3 dist3 = Distance(posRadA, posRadB);
-				real_ rSPH = paramsD.HSML;
 				real_ d = length(dist3);
 				if (d > RESOLUTION_LENGTH_MULT * paramsD.HSML) continue;
+
 				real4 velMasB = FETCH(sortedVelMas, j);
 				real4 rhoPresMuB = FETCH(sortedRhoPreMu, j);
 

__device__ inline real3 DifVelocity_SSI_DEM(
				const real3 & dist3,
				const real_ & d,
				const real_ & rSPH,
				const real4 & velMasA,
				const real4 & velMasB) {
//printf("** DifVelocity_SSI_DEM\n");
	real_ l = paramsD.MULT_INITSPACE * rSPH - d; //penetration distance
	if (l < 0) {
		return R3(0);
	}

//	real_ c_ab = 10 * paramsD.v_Max; //Ma = .1;//sqrt(7.0f * 10000 / ((rhoPresMuA.x + rhoPresMuB.x) / 2.0f));
//	//real_ h = paramsD.HSML;
//	real_ rho = .5f * (rhoPresMuA.x + rhoPresMuB.x);
//	real_ nu = alpha * paramsD.HSML * c_ab / rho;

//	//*** Artificial viscosity type 1.2
//	real_ nu = 22.8f * paramsD.mu0 / 2.0f / (rhoPresMuA.x * rhoPresMuB.x);
//	real3 derivV = -velMasB.w * (
//		rhoPresMuA.y / (rhoPresMuA.x * rhoPresMuA.x) + rhoPresMuB.y / (rhoPresMuB.x * rhoPresMuB.x)
//		- nu * vAB_Dot_rAB / ( d * d + epsilonMutualDistance * paramsD.HSML * paramsD.HSML )
//		) * gradW;
//	return R4(derivV,
//		rhoPresMuA.x * velMasB.w / rhoPresMuB.x * dot(vel_XSPH_A - vel_XSPH_B, gradW));

	//*** Artificial viscosity type 2
	real_ rAB_Dot_GradW = dot(dist3, gradW);
	real_ rAB_Dot_GradW_OverDist = rAB_Dot_GradW / (d * d + epsilonMutualDistance * paramsD.HSML * paramsD.HSML);
	real3 derivV = -velMasB.w * (rhoPresMuA.y / (rhoPresMuA.x * rhoPresMuA.x) + rhoPresMuB.y / (rhoPresMuB.x * rhoPresMuB.x)) * gradW
			+ velMasB.w * (8.0f * multViscosity) * paramsD.mu0 * pow(rhoPresMuA.x + rhoPresMuB.x, -2) * rAB_Dot_GradW_OverDist
					* R3(velMasA - velMasB);
