Comments for post 2010-12-15-why-my-fluids-dont-flow

[tommadams]

Comments for blog post 2010-12-15-why-my-fluids-dont-flow

[tommadams]

comment imported from WordPress
Niall said on 2011-01-05 21:17:48:

Hey, I'm the guy that just left a comment on one of your videos.
Excuse me if I go on a bit of a ramble, you've been experiencing the same sort of issues & discoveries I have recently!

I fully implemented the viscoelastics paper a few weeks ago in flash - http://niallmeister.deviantart.com/art/Viscoelastic-Fluid-Sim-174243339 - albeit with a few small changes. I've since redone it quite a bit, switched to Leapfrog integration and completely changed the spring algorithm to suit (hey, that whole bit isn't even slightly physically accurate anyway, as long as it's stable right?).

The particle clumping effect you're describing is pretty well known, referred to as a tensile instability. There have been various efforts to reduce / eliminate the effect, a good paper on it is "SPH Without a Tensile Instability" by J.J. Monaghan. I agree though, the results you get from double density relaxation are far nicer looking and more stable in general.

There's a common, very annoying issue specifically with the "Müller-esque" formulation of SPH. In 2 fluids of a high density contrast a gap tends to form between them, it's detailed and a fix is proposed in the paper Density Contrast SPH Interfaces http://www.ifi.uzh.ch/arvo/vmml/admin/upload/Solenthaler_sca08.pdf
I've tried doing multiphase sims with double density relaxation, things like changing the pressure force to massi*massj*(pressi+pressj)*weight+(pressneari+pressnearj)*weight*weight and so on but the gap starts to appear there too (and it's not very physically viable, I don't particularly mind though :P). Your implementation is pretty clever, I'm just wondering how you arrived at the kNorm and kNearNorm variables?

[tommadams]

comment imported from WordPress
Niall said on 2011-01-05 21:43:35:

Sorry, wrong link, that's an even older version, here's the correct one http://niallmeister.deviantart.com/art/Viscoelastic-Jelly-183686746

Also sorry for inadvertently plugging my stuff, haha.

[tommadams]

comment imported from WordPress
Tom Madams said on 2011-01-05 22:45:47:

Hi Niall, it's great to hear from people working on the same stuff. I'd definitely noticed both the problems you mentioned, but hadn't come across either of those papers, to thanks for the links. I'm looking forward to checking them out.

Congrats on the Flash simulations, I'm impressed that the performance is so good. You've motivated me to try and implement the full viscoelastic sim.

The kNear and kNearNorm variables are the normalization constants for the two kernel functions. They're the scaling values required to ensure that the kernels integrate to 1.

If I remember correctly, the normalization constants are folded directly into the stiffness and rest pressure constants in the viscoelastic paper. I guess that's in keeping with the tone of that paper (keeping things simple), but I felt like it was hiding important details about the SPH method in general, so reinserted them explicitly back into my code.

When analytically integrating the kernel functions in 2D, you have to add an extra 2*pi*r term to account for the fact that the functions are specified in a circular coordinate system.

Given the kernel function:
(1-r/h)^3

Integrating from 0 to h in circular coordinates becomes:
Integrate[((1-r/h)^3) * r*2*pi, {r, 0, h}]
Here's the result courtesy of wolframalpha.com

Which evaluates to:
(2*pi*h^2)/20

Dividing by this constant ensures the kernel is normalized.

kNearNorm is similarly derived.

Since most of the papers talk about SPH in 3D, they're working in a spherical coordinate system, so the integration is different. That caught me out for the longest time.

[tommadams]

comment imported from WordPress
Niall said on 2011-01-06 02:23:28:

Thanks for the speedy elaboration! I'm only 13 so the maths side of this stuff is still pretty new to me and I haven't fully dived into integrals yet, I somewhat understand though.
One thing I'm wondering about, it doesn't seem like when integrating the equations of motion for the particles you don't take into account their mass, for acceleration. Am I overlooking it or has it been omitted?

[tommadams]

comment imported from WordPress
Tom Madams said on 2011-01-06 11:04:35:

You're right, I apparently forgot to account for mass in the viscosity force calculation.

As for surface tension and pressure relaxation, I actually calculate the acceleration directly (rather than the force), which is where the division by pi.m comes from in lines 275, 276, 281 and 282.

That being said, it's entirely possible that I've made a mistake somewhere as my maths isn't terribly good either :/

By the way, there were a couple of mistakes in my previous reply, which I've now corrected. Most notably, the extra integration term should be 2*pi*r, not pi*r.

[tommadams]

comment imported from WordPress
Niall said on 2011-01-08 12:40:41:

Hey, I've implemented your method of doing things myself in flash. I haven't experimented with it much as I'm posting this comment about 10 minutes after getting a working build because I've noticed an issue. The mass/density contrast, aargh!

Here's a demo with a timestep of .1, with 2 fluids of mass 1 and 30 respectively - http://www.fileize.com/view/c82523d7-915/ - mouse to interact, sorry that it's a bit like watching paint dry at the moment.
I also have yet to try varying rest densities, that could possibly contribute to the effect as well.

Do you think I might've coded this wrong or have you experienced the same thing?

[tommadams]

comment imported from WordPress
kaoD said on 2011-02-08 11:03:20:

Beautiful videos! I got all scared of SPH when I faced it the first time, I might give it a try.

I tried your code to see if I could make a video generator out of it, for the sake of entertainmient, but particle radius below 0.04f makes it crash. How do you make the videos then, if it crashes for small radius? With a huge rendering area?

[tommadams]

comment imported from WordPress
Tom Madams said on 2011-02-08 11:13:57:

I never got any hard crashes when running it locally. Then again, the code I posted is a much simplified version of what I actually use at home, so it's possible I introduced a bug somewhere.

In the posted videos, the particle radius is unchanged at 0.05 and kViewWidth is increased to something like 80 or more if I remember correctly, which zooms the camera out.

I hope this helps; I'd love to see your videos if you come up with anything cool.

[tommadams]

comment imported from WordPress
kaoD said on 2011-02-08 13:13:32:

With crash I mean unconsistent behaviour. All particles explode and go crazy.

Anyways I think you might've introduced a bug, everything tends to splash a lot, even sudden splashes out of nowhere at the bottom of the field. Particles tend to flow like marbles or sand too. I'll play with constants and see if I can fix it.

I guess the part you simplified is the drawing part, right? I can clearly see the spaces between adjacent particles, so I guess the real code does some smoothing.

[tommadams]

comment imported from WordPress
Tom Madams said on 2011-02-08 13:36:08:

Ah that's not a bug, it's what I meant by "Unfortunately, SPH simulations have a tendency to explode if the time step is too large" - it's caused by forces in the simulation being too large for the integration method.

You could replace the simple forward Euler with some more complicated integration scheme, but the simplest thing to do is to increase the number of substeps per frame. For some of the videos, I was using around 20 substeps.

Alternatively, if you don't mind your fluid being a little uh, compressible, you could also reduce the stiffness constants.

[tommadams]

comment imported from WordPress
kaoD said on 2011-02-12 08:11:12:

Nice, got it to work and had some fun tweaking constants, but it takes so long per frame after a large number of particles is spawned, and re-running it from the beggining after each tweak is a bit time consuming.

I'm planning on doing this realtime or at least a bit faster, implementing a quadtree for the neighbor check and Leapfrog integration to reduce substeps. If I succeed in doing so I will let you know :)

By the way, was I right when I guessed that the simplified part is the drawing one? A big number of particles would make it look continous like in your videos, or it had some kind of near-deformation to make it look liquid?

[tommadams]

comment imported from WordPress
Tom Madams said on 2011-02-12 10:46:59:

I'm looking forward to seeing what you come up with :)

The rendering code is actually more or less unchanged - you could always increase the size of the constant in glPointSize if you feel the points are too small.

The major simplification in the code I posted is the removal of multi-threading.

The code as posted is set up to run on multiple threads; if you're only running on one thread, you could rearrange it and probably speed it up by around 30%. For example, applying forces to both particles involved in an interaction at the same time.

[tommadams]

comment imported from WordPress
Mikhail said on 2011-06-21 08:21:05:

Hello, Tom!
I would like to talk with you about SPH fluids, could you drop me your skype?

[tommadams]

comment imported from WordPress
Tom Madams said on 2011-06-21 16:53:16:

Hi Mikhall, I'm afraid I don't have Skype (or much free time these days). I've dropped you a mail.

[tommadams]

comment imported from WordPress
Mikhail said on 2011-06-24 08:08:43:

Thanks! I was sent some questions to you via e-mail, hope you'll react as soon as you can.

[tommadams]

comment imported from WordPress
yanioaioan said on 2012-04-03 09:10:31:

Hi there,

i was lookin gat your code above and i can't seem to find where do u implement the Tait's equation of state to calculate pressure (as it was done in 2007 paper by Becker et. al ). I am having an sph algorithm where i firstly calculate my densities based on neighbour search and then i try to use Tait's EOS so as to calculate pressure but it's explodes. i see you use the simple gas law as i did.

Any ideas to share would be great!

Cheers

[tommadams]

comment imported from WordPress
Tom Madams said on 2012-04-04 18:47:14:

My code is based mostly on the Particle-based Viscoelastic Fluid Simulation paper by Clavet et al, which IIRC doesn't use Tait's equation. Instead they opted for more a approximate model that tends to explode less.

[tommadams]

comment imported from WordPress
Jumes said on 2013-06-05 16:07:14:

Hi,
I run your code, then I find there is large space between each particle.It seems that it is just a circle.Your demo video looks better.Building surface tringle mesh? How did you render the water surface?

[tommadams]

comment imported from WordPress
Tom Madams said on 2013-06-06 08:21:29:

For the most part, the videos are rendered using many more particles than in the demo.

To improve upon simple circles for rendering the particles, one common technique is to render gaussian splats for each particle into a off-screen accumulation buffer and then perform a per-pixel thresholding operation when rendering the accumulation buffer to the back buffer. Another technique that was used PixelJunk Shooter to hide some of the large gaps (especially when the fluid was falling) was to decay the accumulation buffer by some amount each frame, rather than clearing it.

I don't have any references to hand for either of these techniques, but hopefully that should give you pointers.

[tommadams]

comment imported from WordPress
Bakatare said on 2013-07-07 05:23:49:

This may sound awfully retarded, but what exactly is the gradient of the kernel function in the very first equation here (e. g. nabla W)? I mean, W is not a scalar field, just a simple function of one argument (the distance). Must be some gap in my math, have read several papers on SPH, it's used everywhere and I totally don't get it. Could you, please, explain how to compute that?

[tommadams]

comment imported from WordPress
Tom Madams said on 2013-07-08 10:42:13:

Oh dear, I'm afraid it's been a bit too long since I wrote this post and I've forgotten all the details :(
I didn't actually use the pressure model that the first equation describes, (I use the model from Particle-based Viscoelastic Fluid Simulation) so that gradient doesn't appear in the source code I provided, which is why that might be confusing.

The the thing to remember is that although W might be defined in terms of distance from particle i, it still forms a scalar field. Think about a particle i in 2D space: W(r_i, r_j, h) forms a non-zero scalar field in the region of space described by a circle of radius h centered at r_i. I hope this helps at least a little.

[tommadams]

comment imported from WordPress
Me said on 2013-09-17 20:13:13:

Hi Tom,

I just wanted to say thanks for sharing your code. It helped me a lot getting my own SPH working.

http://softologyblog.wordpress.com/2013/09/17/multiphase-smoothed-particle-hydrodynamics/

Regards,
Jason

[tommadams]

comment imported from WordPress
Tom Madams said on 2013-09-17 21:13:23:

No worries, I'm glad you found it some help.

[tommadams]

comment imported from WordPress
Rhody Lugo said on 2014-01-10 03:20:54:

Pretty cool stuff!
I just wanted to ask it your simulation models the incompressiblility of the liquids.
Regards.

[tommadams]

comment imported from WordPress
Tom Madams said on 2014-01-10 09:10:40:

No, the simulation makes several approximations in order to provide more stability at larger time steps. This is why you can see some compressibility in the videos.

[tommadams]

comment imported from WordPress
JPK said on 2014-03-14 13:31:34:

recently im trying to develop my own version of this in python. somehow this code doesnt make use of the standard formulas used in most of the papers. i have a few questions - maybe someone who is interested in sph can contact me so we can solve it together. probably the author of this isnt working on sph for a few years now :D

1. The Kernels - why do they need to be normalized. and why are there often normalization factors for 2d and 3d given? isnt the Kernel basically always a 1d function? some even say they need to be symmetric W(r,h)=W(-r,h) - but as r is the Euclidian distance this doesnt make any sense to me.

2.How can you convert pseudo densities and pressures to real-life quantities? I wonder how you can compute a real density by adding weighted masses together??

3. Someone above postet a similar question regarding the gradient. as i only have a simple function - is it correct that i just calculated the direction vector between particle i and j. and then i multiplied the vector with the result of the W-function (1. derivative). so i have the direction pointing away from the calculated particle and the length is the value of the function. later i added all vectors together to get the final gradient?

4. are there papers regarding compressible fluids like air? ive only seen water :(

[tommadams]

comment imported from WordPress
crococode said on 2014-04-11 07:48:50:

Hey Tom! Thanks a lot for your code posted here. I used your basic layout to get my implementation started. Still fighting with some surface tension issues for varying particle sizes, but the first results were really good:

http://crococode.wordpress.com/2014/04/11/smoothed-particle-hydrodynamics-even-more/

[tommadams]

comment imported from WordPress
hassan said on 2014-04-11 08:04:39:

Hi
Have you used the kernel gradient correction term in your code? please

[tommadams]

comment imported from WordPress
Golzar said on 2014-06-05 07:43:15:

Hi Tom
For interaction between two phase (different particles with different density ) what did you do it?
Regards