Integration of FLUIDS v.2 (SPH Fluids)

[Rvo]

Apologies, seems I've made a mistake with copying the new cpp files for the solvers or something as it now works indeed.
I couldn't add particles to the fluid as it would crash and some weird stuff...

One thing though, when I want to use the GPU solver and get the particle position from the fluid and substract a contact point location, I get 1.#QNAN0 as vector values (99% anyway). How does one fix that to get floats/doubles or am I doing something wrong?

Code: Select all

btVector3 particlePosition = getFluid()->getPosition(m_particles[particleIndex]->m_fluidParticleIndex) - m_particles[particleIndex]->m_lastHitPoint

Seems to work with the CPU solver so hmm.

[rtrius]

The OpenCL solver only generates a force and does not directly modify the positions of
the particles(aside from sorting values during the grid update), so the most likely cause
is a NaN SPH force, which would be propagated to the velocity and then position.

Is the velocity also NaN?

Try looking at
FluidSphDemo\BulletMultiThreaded\SphSolverOpenCL\btFluidSphSolverOpenCL.cpp and
FluidSphDemo\BulletMultiThreaded\SphSolverOpenCL\btFluidSphOpenCL.h

In order to determine whether the solver is giving a NaN SPH force, it would be worth
looking at 2 points in btFluidSphSolverOpenCL.cpp.

The first is near line 163. In order to access the position, velocity, etc. it is necessary to
first transfer the data from GPU to CPU. btOpenCLArray is used to store data on the GPU,
and btOpenCLArray::copyToHost() can be used to copy data back to the CPU.

For example:

Code: Select all

btFluidSphOpenCL* fluidData = ...;

static btAlignedObjectArray<btVector3> position;
fluidData->m_pos.copyToHost(position, true);

The second is at line 217; m_tempSphForce contains the force that has been read back to the CPU.

If there are NaNs near either of these points then the issue is probably caused by the solver.

It might also be worth setting UPDATE_GRID_ON_GPU on line 95 to false, and checking whether
that resolves the issue.

[Rvo]

Alright, I'll try that, thanks.

One more thing.. is it possible to ask BulletFluids how much particles are with a certain radius of a particular particle? When a particle is hanging from my mesh, I want to kinda see if there are 4-5 particles close to it, so they would "clump" together and scale down the amount of force I apply to make the particles stick to the mesh. So its kind of a fake version of multiple particles together forming a larger drop that is heavier and might be to heavy to keep it sticking. So its basically an int returning function saying how many particles are in the vicinity (radius A) of particle B. I'm not sure this is actually easy to do.. but if you might have a good suggestion that would be very welcome

Trying to get the closest point on my mesh in respect to the particle's position but that seems a bit harder than I thought. I figured maybe I could ask a solver to find the nearest point on my collision shape but not much luck. The GJK-algorithms in the demos won't work either since one of my shapes is a BvhTriangleShape and the particle? No idea. Raycasting seems like my "best" option but then I'd still have to know which direction to cast my rays, instead of "everywhere". Guess here that expansion radius really comes in handy (:

[dern23]

Rtrius, nice looking work here! Just one question: I'm having problems achieving the framerate levels seen in the youtube link. I've tried running the code on a machine with a GTX 570 and one with a GTX 690, and in both cases the framerates were in the low single digits for the fluids demo. All I did was check the code straight out of Github, used the vs2010.bat file to generate the solution files etc., and set the startup project to App_Bullet3_OpenCL_Demos_NVIDIA. I'm compiling in release mode with VS2010. I know there's got to be something going on, but I'm at a loss! Any suggestions would be greatly appreciated

[rtrius]

One more thing.. is it possible to ask BulletFluids how much particles are with a certain radius of a particular particle?

If using the default or multithreaded CPU solver, the neighbor tables contain indices and distances
of particles within the SPH interaction radius(typically 20-40 particles). The number of particles is
not explicitly stored so a loop would be needed to get it. Additionally, the distance is at simulation
scale, so make sure to convert it to world scale before using it.

To access the neighbor table(see btFluidSphSolver.h for struct defines):

Code: Select all

btFluidSphSolverDefault::SphParticles* solverData = static_cast<btFluidSphSolverDefault::SphParticles*>( btFluidSph::getSolverData() );

To get the number of neighbors:

Code: Select all

btAlignedObjectArray<int> numParticlesInSphRadius;	//1 entry per particle

for(each particle)	//index i
{
	const btFluidSphNeighbors& table = ...;
	for(each entry in neighbor table)
	{
		int n = table.getNeighborIndex(...);

		//Each entry in the neighbor table contains a pair of particles
		//(for neighbor table of particle i, index i is implicit)
		numParticlesInSphRadius[i] += 1;
		numParticlesInSphRadius[n] += 1;
	}
}

In order to get the number of particles within a radius greater than the SPH interaction radius, it would
be necessary to create another grid and recompute distances; much slower.

The OpenCL solver has no built in feature for this, so modifications would be needed(such as adding an
additional kernel) or manually recalculating the distance on CPU as above. It also does not really begin
to outperform the multithreaded solver(4 threads), until ~16000+ particles are used.

Just one question: I'm having problems achieving the framerate levels seen in the youtube link. I've tried running the code on a machine with a GTX 570 and one with a GTX 690, and in both cases the framerates were in the low single digits for the fluids demo.

I'm not sure why the performance is so low.

Could you pull the latest changes(just disabled some printf), rerun premake, run the demo in profiling
mode and attach the error.log here?

From the Command Prompt, navigate to \bullet3\bin\, and type:

Code: Select all

\bullet3\bin\App_Bullet3_OpenCL_Demos_NVIDIA_vs2010.exe --x_dim=0 --y_dim=0 --z_dim=0 --selected_demo=12 --dump_timings >> error.log

 --x/y/z_dim = 0 to start with 0 rigid bodies,
 --selected_demo = 12 to start with the fluid demo, and
 --dump_timings to output profiling information.

[Rvo]

Thanks for the neighboring particle example, I will try to implement it!

In the manual pdf that comes with the library, you mention the following:

Note that, as this is a simulation scale parameter, the default gravity is 125x larger compared to rigid bodies. This default value can be reduced considerably. (

What's a good value here? I'm setting it to -9.8 for the fluid world and don't change anything for the local parameters. I think the default there is -4.8. Will it act realistic like that?
------------------------

I'm updating my contacts every internal callback. You set the distance for each contact, in the code it is "distance + m_fluidLocalParameters.m_particleRadiusExpansion". This distance quite differs from my own calculation where I do the following:

Code: Select all

btVector3 distanceVector = getFluid()->getPosition(m_particles[particleIndex]->m_fluidParticleIndex) - contact.m_hitPointWorldOnObject;
double realDistance = distanceVector.length();

This gives me a calculated distance of (real example) 0.00422436 (distracted my 0.025f particle radius) versus -0.000224192 according to the contact.m_distance. This kinda confuses me. I use the expansion radius to look around the particle. If it's disconnected from my mesh, I take the X, Y and Z offsets and calculate forces to push it back to the nearest point on the mesh, which would be possible with the expansion radius. Determining the distance to that point could be done with contact.m_distance. But since I need the vector components, I split up distanceVector in X, Y and Z, like in my code example. The order I do it in, in the post tick internal callback: update indices for particles, update contacts, update visual particle positions and then apply my blood/sticky forces things.

Doing the last hitpoint calculation manually (posttick , while getting latest contacts) gets me much better distances, very very close to zero:

Code: Select all

_physicsEngine->m_particles[i]->m_lastHitPoint = m_fluidSph->getPosition(contact.m_fluidParticleIndex) - (contact.m_normalOnObject* PARTICLERADIUS);

However at the edges of the mesh when the particle flies off, this hitpoint is actually not precisely on the edge... which seems all very strange to me. Using my own code I get better(?) up to date hitpoints that actually give ~0.00001 distance from particle to mesh instead of 0.004-ish values. Although I need to account for expansion radius still... what do you make of this, if you understand what I typed up here (sorry for the wall of text).

If I take the contact's last hitpoint (in blue) then this is what happens each frame: particle is in front of the contact somehow. My code does kinda "fix" that, somehow (second image).
But yeah as you can see, the edge detection is a bit off, being the last hitpoint before the particle flew off?



I guess I kinda need the real contact point, without any expansion to know how far the particle is from the mesh. Then, if the distance from the particle to the mesh is too great (particle position vs last known hitpoint), consult the expansion radius to find the nearest point in the mesh. Otherwise I don't know how to tackle this exactly. I cant go raycasting every frame into X, Y and Z to see if I encounter my mesh there... Any ideas?

------------------------
I'm trying to create bleeding wounds with your library. I've applied a normal distribution to the spawning of particles at the emitter. So the mean for the distribution is halfway between 0 and the particle max_lifetime, after which it is relocated to the emitter. Anyway. At one point in time I will have "many" particles (10-20) spawning very close to each other. So they will actually very likely overlap each other. Now with the stiffness at 0.5 they will continue to push each other away which looks pretty bad. Yet it is blood so the stiffness should be somewhere around 1.5 actually, being in-compressible and all? Do you have any suggestions or should I regulate the spawning per emitter so that the particles will not overlap (and push each other away) and then maybe place several emitters close to each other to simulate a cut or bullet hole? I kinda hoped to make 1 emitter one bullet hole but that proves a bit hard to accomplish.

I guess another option is to spawn particles not AT the emitter but around it, with an offset matching the particle radius. Then I'd have to look at other emitters though, making it impossible to place then nearer than 4* particle radius to each other (making room for 2 simultaneously spawning particles at 2 different emitters).
------------------------


------------------------
For the sake of my experiment, I'd like the particles to maybe look more like they can group up and form a larger droplet. However, there's a spacing between particles now that I'm not totally sure if that can be removed. I do this:

Code: Select all

FG.m_simulationScale *= (1/PARTICLE_RADIUS);

My particle radius is 0.025f, so the simulation scale here 1/(0.004*40) = 6.25. That works fine. Now, I think I know that the fluid should be rendered via some different technique than just creating a mesh/billboard at a particle position. However, should this be the desired rendering technique, is it possible to move the particles closer to eachother? My billboards/meshes are 2* particle radius but they aren't actually touching at the moment, which I'd like for visual purposes. The grid cellsize is 0.0625. I'm pretty sure the stiffness doesn't really matter in this case? Am I overlooking something or is this actually not an intentional use of the library? I kinda doubt the latter.



I can ofcourse make the mesh/billboard 2* larger so they will kinda look more attached but that's kinda besides the point maybe?

[dern23]

Hmm, when I ran just the fluid with no rigid bodies it was very fast. I also got the same slowdown in the other demos as well on both of the GPUs I listed, so the problem is probably Bullet related rather than anything specific to your code. I'll just make another topic to avoid hijacking this one. Thanks anyway!

[rtrius]

What's a good value here? I'm setting it to -9.8 for the fluid world and don't change anything for the local parameters.
I think the default there is -4.8. Will it act realistic like that?

-9.8 gravity is fine; whether it is realistic depends on the definition of 'realistic'. With the default
parameters, a gravity of -1.9 is needed to make a column of fluid to fall at the same rate as a
rigid body( world scale gravity -9.8 ). At the same time, though, that gravity also makes the fluid
much more floaty, as if the simulation were performed on the Moon instead of the Earth.

When performing a real time simulation, especially on CPU, there are not really any 'more realistic'
or 'more correct' values due to the amount of approximations/hacks. For example, an accurate
stiffness for water is on the order of 10^6 or 10^7, which would require time steps less than 10^-6.
The right way to simulate an incompressible fluid with SPH is to use an iterative method (IISPH, PBF,
or constraint fluids), but those methods are also more expensive so fewer particles can be used. On
the CPU, it would amount to using, say, 3000 particles instead of 16000 particles.

I'm updating my contacts every internal callback. You set the distance for each contact,
in the code it is "distance + m_fluidLocalParameters.m_particleRadiusExpansion". This distance quite differs from my own calculation where I do the following...

The expansion radius is already accounted for; that is, the distance in the contact does not include the
expansion radius.

If there is a particle with radius 1 and expansion radius 0.5, then a sphere of radius 1.5 is used for
the collision detection. Say that there is a plane at y = 0.0, and a particle as above at y = 1.5. The
collision detection returns a distance of 0.0, but the actual distance is 0.5. For this reason, the
distance is calculated as "distance + m_fluidLocalParameters.m_particleRadiusExpansion".

The contact point being a bit off is due to the position update; the internal loop is:

• -Calculate and apply SPH forces
  -Perform collision detection(generate contacts)
  -Correct velocity
  -Update position

So the contacts(hit point, distance, and normal) are off by 1 frame. A single frame should not make
much of a difference if a force is applied to move the particles towards the mesh, but
btFluidSphRigidCollisionDetector::performNarrowphase() can be used to update the contacts if the
current state is absolutely needed. (If that is too slow, then it would be necessary to add another
callback.)

Anyway. At one point in time I will have "many" particles (10-20) spawning very close to each other. So they will actually very likely overlap each other. Now with the stiffness at 0.5 they will continue to push each other away which looks pretty bad. Yet it is blood so the stiffness should be somewhere around 1.5 actually, being in-compressible and all? Do you have any suggestions or should I regulate the spawning per emitter so that the particles will not overlap (and push each other away) and then maybe place several emitters close to each other to simulate a cut or bullet hole?

Not sure about the emitters, but it would probably help to use a lower stiffness and m_initialSum to
make the particles push on each other less. The default values(1.5 and 0.25) are both very high.
0.5 for the stiffness should be fine, try reducing the m_initialSum.

For the sake of my experiment, I'd like the particles to maybe look more like they can group up and form a larger droplet.
However, there's a spacing between particles now that I'm not totally sure if that can be removed. ...

Surface tension forces are necessary to make the particles form into clumps, but they are not
currently implemented. If I can find the time I'll look into implementing it.

The SPH pressure force tries to maintain a constant density, so to make the particles move closer
you could decrease the SPH particle mass, but it could mean that other parameters need to be
changed. The stiffness actually does have an effect on the particle spacing, but it is not as noticeable
unless there is a column or pool of fluid.

I don't quite understand your last point though; what is the issue with increasing the billboard radius?

[Rvo]

-9.8 gravity is fine; whether it is realistic depends on the definition of 'realistic'. With the default
parameters, a gravity of -1.9 is needed to make a column of fluid to fall at the same rate as a
rigid body( world scale gravity -9.8 ). At the same time, though, that gravity also makes the fluid
much more floaty, as if the simulation were performed on the Moon instead of the Earth. [...]

Fair enough, the "floatyness" was exactly my point. With -9.8 it seems that particles are pulled downwards too hard (although it actually might be correct, droplets do tend to fall pretty fast), with -4.8 they float to the ground. I prefer the -9.8 value I think.

The contact point being a bit off is due to the position update; the internal loop is:

• -Calculate and apply SPH forces
  -Perform collision detection(generate contacts)
  -Correct velocity
  -Update position

So the contacts(hit point, distance, and normal) are off by 1 frame. A single frame should not make
much of a difference if a force is applied to move the particles towards the mesh, but
btFluidSphRigidCollisionDetector::performNarrowphase() can be used to update the contacts if the
current state is absolutely needed. (If that is too slow, then it would be necessary to add another
callback.)

Not sure about this one. Doing some tests it feels like its much more precise to have the exact contacts at each frame. When a particle is rolling over the mesh surface, the distance between the contact and the particle is around 0.004(m). When the particle disconnects from the mesh, I want to detect minimal distances so that I can correct them right away. Therefore I'm thinking that the distance between the particle and the contact are always off by 0.004m (in my current setup) - which isn't exactly much but it's already ~1/6th of my particle radius.

I'm not totally sure what you mean by using btFluidSphRigidCollisionDetector::performNarrowphase() - do you mean a call to it? I cant seem to access m_fluidRigidCollisionDetector.performNarrowphase() directy. Or do you mean adding code to re-evaluate the contacts? Same for the callback, do you mean an internal pre/post tick callback or? Sorry for the misunderstanding.

I don't quite understand your last point though; what is the issue with increasing the billboard radius?

Well, the billboard radius matches the particle radius. I can increase the radius but then the particles do not match the billboard anymore. So that would mean billboards intersecting with the mesh, which kinda looks bad. So that's why I'd prefer particles not pushing eachother away so much. I will try your sum suggestion (: I know surface tension etc. would be a lot more realistic but I'm looking for a low-cost 'solution', in this case the sticky force in the Y-direction is scaled down when more neighbors are detected. This would "simulate" the accumulation of particles, creating larger drops, making it heavier.

Any ideas for using the expansion radius? I'd expect when a particle is close to the mesh, rolling over it even, that the contact is below the particle, even if the expansion radius is quite large (0.5f)? In my current setup I get that but after a while the new contact will be ahead of the particle. It makes sense, but its not exactly working for me atm:

What I'd really want, for my purposes that is, is keeping the contact under the particle until the edge, then the particle will fly off and the expansion radius will now now detect the mesh still as the particle falls down the side. I hope you know what I mean. Is there any solution to this? If not, do you know how I can ask a particle to determine its closest rigid body that I want it to stick to? I can't get it to work with regular Bullet methods, so you're kind of my only hope here. I'm thinking along the lines of using the distance without expansion until that distance is > some value, probably the particle radius (so not touching mesh anymore) and then switching to the distance using the expansion...? Or if the particle has contact but suddenly stops and use the contact.m_distance? I'm not sure.

Thanks!

[rtrius]

I'm not totally sure what you mean by using btFluidSphRigidCollisionDetector::performNarrowphase() - do you mean a call to it? I cant seem to access m_fluidRigidCollisionDetector.performNarrowphase() directy. Or do you mean adding code to re-evaluate the contacts? Same for the callback, do you mean an internal pre/post tick callback or? Sorry for the misunderstanding.

No need to apologize, I probably could have written it more clearly. Creating a new
btFluidSphRigidCollisionDetector, and then calling performNarrowphase() (with the btFluidRigidWorld's
dispatcher, and dispatchInfo) can be used to re-evaluate the contacts with the current particle positions.
The btFluidSphRigidCollisionDetector does not store any data or state, so creating another should
not cause any side effects.

As for the callback, I mean to add a new 'mid-tick' callback in between the pre and post tick callbacks
that is called right when the forces are applied and the contacts are current.

Any ideas for using the expansion radius? I'd expect when a particle is close to the mesh, rolling over it even, that the contact is below the particle, even if the expansion radius is quite large (0.5f)? In my current setup I get that but after a while the new contact will be ahead of the particle. It makes sense, but its not exactly working for me atm:

The reason that the contact is ahead is likely due to the sphere-triangle collision algorithm:

• -Project sphere center onto triangle plane
  -Check if the projected point in the plane is inside the triangle
  -If sphere center is inside, use point-plane distance
  -If sphere center is outside use point-line segment distance(use the nearest distance from the 3 triangle edges)

I think whats happening in that picture is that the expanded radius sphere is colliding with
the edge of the triangle. However, there should be 1 contact for each colliding particle-triangle
pair. If the expanded radius sphere is intersecting with 2 triangles in that picture, then there should
be 2 contacts.

In the current implementation, the contacts in each btFluidSphRigidContactGroup should also
be sorted by the fluid particle index, although it is not explicitly guaranteed by the API.
If a single btRigidBody/btCollisionObject is used for your mesh, it should be possible to
loop through the contact group and, for each particle, find the contact with the closest distance
to get the contact 'below' the particle.

[Rvo]

I feel a bit like a tool. Ofcourse a particle has more contacts than 1, potentially. I'm now taking the closest contact and it works wonderfully... Thanks for bringing that up...

[Rvo]

Hm. In the PostTick I now do this, right before I loop through the RigidContacts to match particles:

Code: Select all

btFluidSphRigidCollisionDetector* newDetector = new btFluidSphRigidCollisionDetector();
const btFluidSphParametersGlobal FG = _physicsEngine->getFluidWorld()->getGlobalParameters();
const btDispatcherInfo dispatchInfo = _physicsEngine->getFluidWorld()->getDispatchInfo();
newDetector->performNarrowphase(_physicsEngine->getFluidWorld()->getDispatcher(), dispatchInfo, FG, m_fluidSph);

However the contact.m_distance is still ~0.005 less than deducting the particle position vector and the hitpoint vector. This is also the case when I do this (0.025 = part.radius):

Code: Select all

(m_fluidSph->getPosition(_physicsEngine->m_particles[i]->m_fluidParticleIndex)-contact.m_hitPointWorldOnObject).length()-0.025

Then the contact.m_distance = 0.005 smaller than the calculated distance. Is this because of penetration? If so, should I avoid using the contact.m_distance and just calculate the diff. between the particle position and the hitpoint and run with that?

[rtrius]

Forgot to mention that btFluidSph::internalClearRigidContacts() should be called before
performNarrowphase(). The updated contacts would be appended at the end of the array
if it is not called, which might be causing the issue.

I'm having a bit of difficulty understanding the question. What is the 'calculated distance'?

Also, the expression:

Code: Select all

distance = (m_fluidSph->getPosition() - contact.m_hitPointWorldOnObject).length() - particleRadius

is exactly how the sphere-triangle distance is calculated, with contact.m_hitPointWorldOnObject
always being the nearest point on the triangle to the sphere center.

[Rvo]

I noticed that too about the array, and I did try and call the clear function although I might have used the non-internal clear (at work atm so can't check) and what I noticed is that I had no contacts left at all when I tried it. Will try it again and check which function I was calling - and if it works you'll not get a peep out of of me.

It is probable I made a mistake with my own distance calculations by just subtracting the particle radius from the (particlePosition - particleHitpoint) vector's length(). I changed this by determining the vector to the particle edge and taking the length() of that vector which seems to be more accurate.

[Rvo]

Doing it like this I get no hitpoints at all. I looked at the narrowPhase code and I didn't see right away why this is the case...

Code: Select all

btFluidSphRigidCollisionDetector* newDetector = new btFluidSphRigidCollisionDetector();
btFluidSphParametersGlobal FG = _physicsEngine->getFluidWorld()->getGlobalParameters();
btDispatcherInfo dispatchInfo = _physicsEngine->getFluidWorld()->getDispatchInfo();
m_fluidSph->internalClearRigidContacts();
newDetector->performNarrowphase(_physicsEngine->getFluidWorld()->getDispatcher(), dispatchInfo, FG, m_fluidSph);

Also, any chance you could work out the neighbor table example some more? I see you have been busy with surface tension but what I'd want for my project right now is an easy fix, which is just summing up the distance-vectors to neighboring particles for each individual particle so I know what average direction the particle is attracted (so particle clumping I guess). I see you can ask distances to neighbors from the table but I kinda need more than that. As for the example, I now get errors already on the first line, "invalid type conversion". I know it can't be difficult but somehow I just can't figure it out...

Code: Select all

btFluidSphSolverDefault::SphParticles* solverData = static_cast<btFluidSphSolverDefault::SphParticles*>(m_fluidSph->getSolverData());

btAlignedObjectArray<int> numParticlesInSphRadius;   //1 entry per particle
for(each particle)   //index i
{
   const btFluidSphNeighbors& table = ...;  <-- btAlignedObjectArray<btFluidSphNeighbors> m_neighborTable; from the Particles struct per... particle?
   for(each entry in neighbor table)
   {
      int n = table.getNeighborIndex(...); <--- ... = i?

      //Each entry in the neighbor table contains a pair of particles
      //(for neighbor table of particle i, index i is implicit)
      numParticlesInSphRadius[i] += 1;
      numParticlesInSphRadius[n] += 1;

      ////////////
      Can I ask distances to neighbors here easily?
      ////////////
   }
}

Thanks and sorry for taking up so much of your time!