 	velMasD[index] = mR4(updatedVelocity, /*rho2 / rhoPresMu.x * */velMas.w); //velMasD updated
 
 	Real4 rhoPresMu = rhoPresMuD[index];
+
+	// 3*** let's tweak a little bit :)
+	if (fabs(derivVelRho.w) > .005 * rhoPresMu.x / dTD) {
+		derivVelRho.w *= (.005 * rhoPresMu.x / dTD) / fabs(derivVelRho.w); //to take care of the sign as well
+	}
+	// 2*** end tweak
 	Real rho2 = rhoPresMu.x + derivVelRho.w * dTD; //rho update. (i.e. rhoPresMu.x), still not wriiten to global matrix
 	rhoPresMu.y = Eos(rho2, rhoPresMu.w);
 	rhoPresMu.x = rho2;

		return;
	}
	Real3 vel_XSPH = vel_XSPH_D[index];
	// 1*** let's tweak a little bit :)
	if (length(vel_XSPH) > .2 * paramsD.HSML / dTD) {
		vel_XSPH *= ( .2 * paramsD.HSML / dTD ) / length(vel_XSPH);
	}
	// 1*** end tweak
	Real3 posRad = posRadD[index];
	Real3 updatedPositon = posRad + vel_XSPH * dTD;
	posRadD[index] = updatedPositon; //posRadD updated

	// *** let's tweak a little bit :)
	Real3 derivV = mR3(derivVelRho);

	derivVelRhoD[originalIndex] = derivVelRho;
}
//--------------------------------------------------------------------------------------------------------------------------------
//calculate particles stresses
