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]