Hi!
Me and some friends have been working on a 2D Smoothed Particle Hydrodynamics implementation in C++ / OpenGL. A few weeks in, we have tried all kinds of approaches, different versions of the equations, different ways to calculate pressure, different ways to set constants and a lot of different kernels. We have tried to implement both viscosity and the XSPH variant, but we always end up with an inplementation where the particles penetrate eachother (leading to high densities and exploding systems).
Here is a short video of our best results so far (1 MB)
http://dl.dropbox.com/u/3639665/sph_wtf2.avi
I was wondering if anyone has any ideas or have had the same problem before. I realize it's hard to go into any deatails beacuse of the amount of code and small things to tweak, but any ideas are much appreciated!
Best regards
Victor
Particle separation in SPH
Judging from the Video, your Particles are barely interacting. Perhaps one of your variables (mass, density, position) is in the wrong dimension?
Could you post a few example Values for what Masses, etc you are using?
Other than that, it could be just a plain coding error. Posting your main SPH functions could help
Could you post a few example Values for what Masses, etc you are using?
Other than that, it could be just a plain coding error. Posting your main SPH functions could help
Author
Our basic setup is as follows:
A Particle class has the data members position, velocity, mass, density, pressure and the density's derivative. It also has an XSPH variable.
We first zero all derivatives, loop through the list, and for every particle we loop trough all the others (no nearest neighbour search yet) where we run the following functions:
There are some constants here and there, I'll get back to how those look like. We made some changes and are trying to tune them in!
kernel, grad and laplacian are constants that we get from constructing the kernels. I'll post those functions later as well!
The pressure does not depend on any of the other particles, so we are updating it after the inner loop:
A Particle class has the data members position, velocity, mass, density, pressure and the density's derivative. It also has an XSPH variable.
We first zero all derivatives, loop through the list, and for every particle we loop trough all the others (no nearest neighbour search yet) where we run the following functions:
There are some constants here and there, I'll get back to how those look like. We made some changes and are trying to tune them in!
kernel, grad and laplacian are constants that we get from constructing the kernels. I'll post those functions later as well!
/* update the density's derivative */
void Particle::updateDensity(const Particle& p,float* grad)
{
dDensity+=p.mass*((velocity[0]-p.velocity[0])*grad[0]+(velocity[1]-p.velocity[1])*grad[1]);
}
/* update acceleration from pressure */
void Particle::updatePressureAcceleration(const Particle& p, float* grad)
{
for(int i=0;i<2;i++)
{
acceleration-=(p.mass*(pressure/(density*density)+p.pressure/(p.density*p.density))*grad);
}
}
/* update acceleration from viscosity */
void Particle::updateViscosityAcceleration(const Particle& p, float laplacian)
{
float sum[2] = {0.0f , 0.0f};
for(int i = 0; i < 2; i++)
{
acceleration += MY*p.mass*((velocity-p.velocity)/(density*p.density))*laplacian;
}
}
/* update XSPH variant */
void Particle::updateXSPH(const Particle& p, float kernel)
{
for(int i = 0; i < 2; i ++)
{
xsph += EPSILON*p.mass*kernel*(p.velocity - velocity)/((density + p.density)/2.0f);
}
}
The pressure does not depend on any of the other particles, so we are updating it after the inner loop:
/* update pressure */
pressure = (PRESSURE_B*(density))*(pow((density/REST_DENSITY),7)-1);
Author
For the kernels, we use this method.
N is the resolution we want in the kernel array, and h is the smoothing length,
Then when we want to use the kernels, we do it like this:
And we use those constants in the functions in my previous post.
We are a bit unsure if we should normalize the kernel or not when we change the smoothing length, we get really different results with different values of h. None are very good, though
N is the resolution we want in the kernel array, and h is the smoothing length,
float* makeCubicSplineKernel(int N, float h)
{
float* kernelConstants = new float[N];
float kernelConst = (10.0/(7.0*PI*h*h));
for(int i=0; i<N; i++)
{
/* q=r/h goes from 0 to 2h, so the last position in the array corresponds to 2h*/
float q = ((float)i*(2.0/N));
if(q <= 1.0) kernelConstants = kernelConst*(1.0-1.5*q*q+0.75*q*q*q);
else if(q <= 2.0) kernelConstants = kernelConst*(0.25*(2.0-q)*(2.0-q)*(2.0-q));
else kernelConstants = 0.0;
}
return kernelConstants;
}
float* makeCubicSplineKernelGradient(int N, float h)
{
float* kernelGradientConstants = new float[N];
float kernelConst = (10/(7*PI*h*h*h));
for(int i = 0; i<N; i ++)
{
float q = ((float)i*(2.0/N));
if(q <= 1.0) kernelGradientConstants = kernelConst*( (9.0/4.0)*q*q - 3.0*q );
else if(q <= 2.0) kernelGradientConstants = kernelConst*( (-3.0/4.0)*(2.0-q)*(2.0-q) );
else kernelGradientConstants = 0.0;
}
return kernelGradientConstants;
}
Then when we want to use the kernels, we do it like this:
// In the main function
float* kernel = makeCubicSplineKernel(KERNEL_N, SMOOTHING_LENGTH);
float* kernelGradient = makeCubicSplineKernelGradient(KERNEL_N, SMOOTHING_LENGTH);
// In the particle neighbour (all the other particles) loop
// i is the current particle and j is the inner loop particle
dx=p.position[0]-p[j].position[0];
dy=p.position[1]-p[j].position[1];
r=sqrtf(dx*dx+dy*dy);
float dxNorm = dx/r;
float dyNorm = dy/r;
// if we are within 2 smoothing lenghts
index = (int)((r/SMOOTHING_LENGTH) * (KERNEL_N / 2.0));
kernelCoeff = kernel[index];
gradientCoeff[0]=(float)kernelGradient[index]*dxNorm;
gradientCoeff[1]=(float)kernelGradient[index]*dyNorm;
And we use those constants in the functions in my previous post.
We are a bit unsure if we should normalize the kernel or not when we change the smoothing length, we get really different results with different values of h. None are very good, though
That cubic spline kernel probably isn't doing you any good. The problem there is that when particles are close together and still moving towards one another, the forces that should be separating them stop increasing. This makes for a fluid that resists compression so weakly (with a high dt) that yeah, it can collapse under its own weight. Especially if the particles all start in nicely-aligned vertical stacks.
Consider using a 'spiky kernel' setup instead. It'll resist that kind of clustering, and it can still be perfectly smooth where 1<q<2, which is where all the nice interparticle interaction takes place anyway.
You can ditch the kernel lookup table, too: sampling a simple spiky kernel uses only basic arithmetic and is always going to be a lot faster than fetching values from a table.
Consider using a 'spiky kernel' setup instead. It'll resist that kind of clustering, and it can still be perfectly smooth where 1<q<2, which is where all the nice interparticle interaction takes place anyway.
You can ditch the kernel lookup table, too: sampling a simple spiky kernel uses only basic arithmetic and is always going to be a lot faster than fetching values from a table.
Author
It's staring to look better in general with three different kernels for density, pressure and viscosity. However, how are we supposed to handle the kernels with regards to the smoothing length? One of the main features of a kernel is that it is supposed to integrate to one. However, none of the defined kernels (Poly6, for example) seems to integrate to one, and also depend on the smoothing length h. Therefore, should we do any kind of normalization with regards to h to be able to change the smoothing length later?
Thanks
Victor
Thanks
Victor
One of the main features of a kernel is that it is supposed to integrate to one.
For an interactive application, that won't matter much at all. It'll make no difference to the user whether the density of the fluid is (on average) 1.0 or 2.0 or 713.25, as you're bound to end up using somewhat arbitrary constants for tuning fluid behaviour anyway, so if your density was too low, you'd just use a higher k until it looked right.
If you do feel you need the kernels to integrate to one, then divide their output values by their integrals. This is just a scalar constant per kernel and only needs to be calculated once - and it can be approximate.
Therefore, should we do any kind of normalization with regards to h to be able to change the smoothing length later?
[/quote]
Nope - just scale the interparticle distance inside the particle loop:
r=sqrtf(dx*dx+dy*dy) / some_smoothing_length_scaling_factor;
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
