Particles Sample

Hi

I was wondering if anyone has developed a sequential version of the Cuda “Particles” sample that ships with the devkit?

I’ve spent some time first coding up a sequential grid version, but found this to be too inefficient. So, I then opted for a kd-tree based solution which seems to work well apart from one aspect that I don’t understand.

I’ve always been fascinated by the “Particles” sample because it doesn’t use complicated rigid-body pairs to handle collisions but yet the results appear correct and impressive.

However, my own sequential mapping of the “Particles” sample follows exactly the same logic and yet I end up with many sphere-sphere interactions with spheres passing right through each other.

There is clearly something I don’t understand in the collision handling and was wondering if anyone else had the same problem or knows what I’m doing wrong.

I’ve never fully understand the particle sorting in the “Particles” sample and was wondering if it is necessary to radially sort the adjacent particles so that nearest collisions are handled first.

Cheers and any help appreciated.

Graham

My own code is pasted below and will be seen similar to “Particles” integrateSystem(), collide() and collideSpheres().

=====

[codebox]

void ParticleSystem::integrate(float deltaTime)

{

for (uint index=0; index<m_numParticles; index++)

{

	float3 pos = mParticles->getParticles()[index].getPosition();

	float3 vel = mParticles->getParticles()[index].getVelocity();

	mParticles->getParticles()[index].setPositionPrevious(pos);

	mParticles->getParticles()[index].setVelocityPrevious(vel);

	float3 gdt = m_params.gravity * deltaTime;

	vel += gdt;

	vel *= m_params.globalDamping;

	// new position = old position + velocity * deltaTime

	pos += vel * deltaTime;

	// set this to zero to disable collisions with cube sides

	if (pos.x > 1.0f - m_params.particleRadius){ pos.x = 1.0f - m_params.particleRadius; vel.x *= m_params.boundaryDamping; }

	if (pos.x < -1.0f + m_params.particleRadius) { pos.x = -1.0f + m_params.particleRadius; vel.x *= m_params.boundaryDamping;}

	if (pos.y > 1.0f - m_params.particleRadius) { pos.y = 1.0f - m_params.particleRadius; vel.y *= m_params.boundaryDamping; }    

	if (pos.y < -1.0f + m_params.particleRadius) { pos.y = -1.0f + m_params.particleRadius; vel.y *= m_params.boundaryDamping;}

	if (pos.z > 1.0f - m_params.particleRadius) { pos.z = 1.0f - m_params.particleRadius; vel.z *= m_params.boundaryDamping; }

	if (pos.z < -1.0f + m_params.particleRadius) { pos.z = -1.0f + m_params.particleRadius; vel.z *= m_params.boundaryDamping;}

	mParticles->getParticles()[index].setPosition(pos);

	mParticles->getParticles()[index].setVelocity(vel);

}

}

void ParticleSystem::collide()

{

// traverse particles

for (uint particleIndex=0; particleIndex<m_numParticles; particleIndex++)

{

	float3 force(0.0f,0.0f,0.0f);

	// collision with adjacent particles from kdtree

	float3 pos    = mParticles->getParticles()[particleIndex].getPosition();

	float3 vel    = mParticles->getParticles()[particleIndex].getVelocity();

	float  particleRadius = mParticles->getParticles()[particleIndex].getRadius();

	float  w2 = 2.0f * particleRadius;

	Point3D<float> lower(pos.x-w2,pos.y-w2,pos.z-w2);

	Point3D<float> upper(pos.x+w2,pos.y+w2,pos.z+w2);

	//int numberAdjParticles = 0;

	for (KDTreeType::iterator kdt_it=mParticlesKDTree->getKDTree()->begin(lower,upper); kdt_it!=mParticlesKDTree->getKDTree()->end(); kdt_it++)

	{

		Point3D<float>	adjPoint    = (*kdt_it).first;

		Particle		adjParticle = (*kdt_it).second;

		int				adjParticleIndex = adjParticle.getIndex();

		if (particleIndex != adjParticleIndex) // dont test against self

		{

			float3 pos2    = mParticles->getParticles()[adjParticleIndex].getPositionPrevious();

			float3 vel2    = mParticles->getParticles()[adjParticleIndex].getVelocityPrevious();

			force += collideSpheres(pos, pos2, vel, vel2, m_params.particleRadius, m_params.particleRadius);

		}

		numberAdjParticles++;

	}

	if (!force.isZero())

	{

		mParticles->getParticles()[particleIndex].updateVelocity(force);

	}

	else

	{

		velCorrections.push_back(float3());

	}

	//std::cout << "numberAdjParticles: " << numberAdjParticles << std::endl;

}

}

// collide two spheres using DEM method

float3 ParticleSystem::collideSpheres(float3 posA, float3 posB,

                                  float3 velA, float3 velB,

								  float radiusA, float radiusB)

{

// calculate relative position

float3 relPos = posB - posA;

float dist2 = lengthSquared(relPos); // just length squared at this stage

float collideDist  = radiusA + radiusB;

float collideDist2 = collideDist * collideDist;

float3 force(0.0f,0.0f,0.0f);

if (dist2 < collideDist2) // particles collide

{

	float dist = std::sqrt(dist2);

    float3 norm = relPos / dist;

	// relative velocity

    float3 relVel = velB - velA;

// relative tangential velocity

    float3 tanVel = relVel - (norm * dot(relVel, norm));

// spring force

	float kd = m_params.spring*(collideDist - dist);

    force = norm * (-kd);

    // dashpot (damping) force

    force += relVel*m_params.damping;

    // tangential shear force

    force += tanVel*m_params.shear;

	// attraction

    force += relPos*m_params.attraction;

}

return force;

}

[/codebox]

Hi

I was wondering if anyone has developed a sequential version of the Cuda “Particles” sample that ships with the devkit?

I’ve spent some time first coding up a sequential grid version, but found this to be too inefficient. So, I then opted for a kd-tree based solution which seems to work well apart from one aspect that I don’t understand.

I’ve always been fascinated by the “Particles” sample because it doesn’t use complicated rigid-body pairs to handle collisions but yet the results appear correct and impressive.

However, my own sequential mapping of the “Particles” sample follows exactly the same logic and yet I end up with many sphere-sphere interactions with spheres passing right through each other.

There is clearly something I don’t understand in the collision handling and was wondering if anyone else had the same problem or knows what I’m doing wrong.

I’ve never fully understand the particle sorting in the “Particles” sample and was wondering if it is necessary to radially sort the adjacent particles so that nearest collisions are handled first.

Cheers and any help appreciated.

Graham

My own code is pasted below and will be seen similar to “Particles” integrateSystem(), collide() and collideSpheres().

=====

[codebox]

void ParticleSystem::integrate(float deltaTime)

{

for (uint index=0; index<m_numParticles; index++)

{

	float3 pos = mParticles->getParticles()[index].getPosition();

	float3 vel = mParticles->getParticles()[index].getVelocity();

	mParticles->getParticles()[index].setPositionPrevious(pos);

	mParticles->getParticles()[index].setVelocityPrevious(vel);

	float3 gdt = m_params.gravity * deltaTime;

	vel += gdt;

	vel *= m_params.globalDamping;

	// new position = old position + velocity * deltaTime

	pos += vel * deltaTime;

	// set this to zero to disable collisions with cube sides

	if (pos.x > 1.0f - m_params.particleRadius){ pos.x = 1.0f - m_params.particleRadius; vel.x *= m_params.boundaryDamping; }

	if (pos.x < -1.0f + m_params.particleRadius) { pos.x = -1.0f + m_params.particleRadius; vel.x *= m_params.boundaryDamping;}

	if (pos.y > 1.0f - m_params.particleRadius) { pos.y = 1.0f - m_params.particleRadius; vel.y *= m_params.boundaryDamping; }    

	if (pos.y < -1.0f + m_params.particleRadius) { pos.y = -1.0f + m_params.particleRadius; vel.y *= m_params.boundaryDamping;}

	if (pos.z > 1.0f - m_params.particleRadius) { pos.z = 1.0f - m_params.particleRadius; vel.z *= m_params.boundaryDamping; }

	if (pos.z < -1.0f + m_params.particleRadius) { pos.z = -1.0f + m_params.particleRadius; vel.z *= m_params.boundaryDamping;}

	mParticles->getParticles()[index].setPosition(pos);

	mParticles->getParticles()[index].setVelocity(vel);

}

}

void ParticleSystem::collide()

{

// traverse particles

for (uint particleIndex=0; particleIndex<m_numParticles; particleIndex++)

{

	float3 force(0.0f,0.0f,0.0f);

	// collision with adjacent particles from kdtree

	float3 pos    = mParticles->getParticles()[particleIndex].getPosition();

	float3 vel    = mParticles->getParticles()[particleIndex].getVelocity();

	float  particleRadius = mParticles->getParticles()[particleIndex].getRadius();

	float  w2 = 2.0f * particleRadius;

	Point3D<float> lower(pos.x-w2,pos.y-w2,pos.z-w2);

	Point3D<float> upper(pos.x+w2,pos.y+w2,pos.z+w2);

	//int numberAdjParticles = 0;

	for (KDTreeType::iterator kdt_it=mParticlesKDTree->getKDTree()->begin(lower,upper); kdt_it!=mParticlesKDTree->getKDTree()->end(); kdt_it++)

	{

		Point3D<float>	adjPoint    = (*kdt_it).first;

		Particle		adjParticle = (*kdt_it).second;

		int				adjParticleIndex = adjParticle.getIndex();

		if (particleIndex != adjParticleIndex) // dont test against self

		{

			float3 pos2    = mParticles->getParticles()[adjParticleIndex].getPositionPrevious();

			float3 vel2    = mParticles->getParticles()[adjParticleIndex].getVelocityPrevious();

			force += collideSpheres(pos, pos2, vel, vel2, m_params.particleRadius, m_params.particleRadius);

		}

		numberAdjParticles++;

	}

	if (!force.isZero())

	{

		mParticles->getParticles()[particleIndex].updateVelocity(force);

	}

	else

	{

		velCorrections.push_back(float3());

	}

	//std::cout << "numberAdjParticles: " << numberAdjParticles << std::endl;

}

}

// collide two spheres using DEM method

float3 ParticleSystem::collideSpheres(float3 posA, float3 posB,

                                  float3 velA, float3 velB,

								  float radiusA, float radiusB)

{

// calculate relative position

float3 relPos = posB - posA;

float dist2 = lengthSquared(relPos); // just length squared at this stage

float collideDist  = radiusA + radiusB;

float collideDist2 = collideDist * collideDist;

float3 force(0.0f,0.0f,0.0f);

if (dist2 < collideDist2) // particles collide

{

	float dist = std::sqrt(dist2);

    float3 norm = relPos / dist;

	// relative velocity

    float3 relVel = velB - velA;

// relative tangential velocity

    float3 tanVel = relVel - (norm * dot(relVel, norm));

// spring force

	float kd = m_params.spring*(collideDist - dist);

    force = norm * (-kd);

    // dashpot (damping) force

    force += relVel*m_params.damping;

    // tangential shear force

    force += tanVel*m_params.shear;

	// attraction

    force += relPos*m_params.attraction;

}

return force;

}

[/codebox]