Water simulation using SPH

Community Forums/Showcase/Water simulation using SPH

Noobody(Posted 2009) [#1]
In the last two weeks I've been working on a water simulation using Smoothed Particle Hydrodynamics.

Implementing this method was quite troublesome, since the derivation of the formulas is a huge chunk of math which was, at least for me, really hard to understand. This is mostly due to most of the maths used not being taught until college. I also had to do a complete rewrite once because I relied too much on a paper on blood simulation using SPH. Said paper used a different set of formulas that are inappropiate for incompressible fluids like water. I then reimplemented the whole algorithm using the approach proposed by Monaghan, which looks quite nice (even though it's more expensive compuationally).

The fluid calculations are more or less optimized and run conveniently fast (the 60 FPS limit is at approximately 7000 particles on my laptop). However, the rendering algorithm is awfully slow. I've already tried to implement it in a shader, but that was even slower. Maybe I'll give it another shot with prerendered metaballs, but I still have to think about that.

Long story short, screenshot:


Edit: Since this thread is getting quite a few hits, I thought I'd post a quick summary of newer and better versions that I posted below this post:
Improved first version
SPH with Soft body coupling
3D version
Videos:
Video 1
Video 2
Video 3

If you still insist on downloading the very first version, here it is:
Download (Win exe + BMax code): Link
/Edit

Fluid particles are added with the left mouse button, boundary particles (=walls) with the right mouse button. Because the rendering algorithm is really slow, the particles are drawn as simple squares by default. You can change the rendering mode with the space bar.

I still have to work on certain discontinuities that express themselves in something that an article calls 'cancer'. For some unknown reason, certain particles suffer from a sudden decrease in density, which induces a huge pressure force that quickly spreads and blows up the whole simulation. I've implemented a hack that finds such particles, constrains them from spreading and tries to convert them back to a regular particle, but it doesn't look right sometimes.

This, this, this and this paper were a great help during the implementation. All formulas I used are derived in those articles.


GW(Posted 2009) [#2]
Pretty darn impressive!


Jeppe Nielsen(Posted 2009) [#3]
Nice work,


plash(Posted 2009) [#4]
Nice implementation.


Shagwana(Posted 2009) [#5]
Very nice, even handles pressure! (U bend style pressure)


BlitzSupport(Posted 2009) [#6]
Yeah, good stuff.


Nate the Great(Posted 2009) [#7]
cool! I personally couldnt find any articles on xsph I dont know how you did but this will help my research considerably. Thanks


markcw(Posted 2009) [#8]
I like.

Could you not optimize it by not calculating particles in the middle of a mass? Or do you do that already?

What's it for?


_Skully(Posted 2009) [#9]
Wow this is fast! I only have dots though... not solid blue like your post


Noobody(Posted 2009) [#10]
Could you not optimize it by not calculating particles in the middle of a mass? Or do you do that already?

Do you mean in the rendering part? I thought about that, but it's hard to determine whether a point is inside or outside the fluid without iterating over all particles and sampling the density at certain points (which the rendering method does - that's what makes it so slow).

What's it for?

Nothing in particular :)
It just fascinated me how realistic the water looked in a few videos on youtube, so I wanted to try if I could make something like this myself.

I only have dots though... not solid blue like your post

See the first post:
Because the rendering algorithm is really slow, the particles are drawn as simple squares by default. You can change the rendering mode with the space bar.



grindalf(Posted 2009) [#11]
Wow thats very impressive


ZJP(Posted 2009) [#12]
Ho yesss. Very nice.
Looking for a blitz3D version now ;) Pleaaaase...

JP


ImaginaryHuman(Posted 2009) [#13]
See my post in the tutorials forum about drawing fast metaballs.


Noobody(Posted 2009) [#14]
I already tried using prerendered metaball textures, but the problem is that those give washed-out borders instead of the sharp boundary I'm looking for. In the second part of the tutorial, you were explaining how to bring out certain treshold values, but those are still a bit washed-out (and also inverted), which doesn't look as good as I want it to be :)

I was thinking about rendering all metaballs in a separate FBO and then drawing that FBO with a glAlphaFunc( GL_GREATER, Treshold ). The problem is that FBOs aren't supported everywhere, so the compability would suffer.

I wonder how programs like OECake work - the water looks extremely good and is also rendered in realtime.


D4NM4N(Posted 2009) [#15]
This is cool! :)


markcw(Posted 2009) [#16]
Do you mean in the rendering part? I thought about that, but it's hard to determine whether a point is inside or outside the fluid without iterating over all particles and sampling the density at certain points (which the rendering method does - that's what makes it so slow).

Yes, in the renderer. Maybe there are iterations that can be done less often?

How about this: determine the largest mass on the screen once every X seconds, then disable particles which aren't moving down any more and are a certain distance from the surface of the largest mass?

Just a thought. If it's not for anything you don't need to optimize anyway.


dawlane(Posted 2009) [#17]
Very nice indeed.


ImaginaryHuman(Posted 2009) [#18]
>>I already tried using prerendered metaball textures, but the problem is that those give washed-out borders instead of the sharp boundary I'm looking for.<<

If you draw rectangles over the whole screen in LIGHTBLEND or SHADBLEND, with a value color of, say, 127,127,127, it will force the outside half of the blob field to be clamped to a threshold. Then do the opposite full-screen rectangle (ie do a lightblend if you did shadeblend, etc), to restore the colors to their original range.

Or something similar.. you could use the stencil buffer to record when you draw your metaballs, then do the above full-screen rects only where the stencil is set, to preserve the background.

You could also clamp the bottom end and the top end of the range so that it leaves all one color value for pixels inside the threshold.


Nate the Great(Posted 2009) [#19]
I wonder how programs like OECake work - the water looks extremely good and is also rendered in realtime.



I believe they use shaders and other things that make them widely incompatible with older machines unfortunately.


DavidDC(Posted 2009) [#20]
Nice, but I get an array bounds crash here
Local Successor:TParticle = FluidGrid[ GridX, GridY ]

when I click near the left window border.

Local GridX:Int = Max(0,Floor( P.PositionX*TSPH.INV_SMOOTHING_LENGTH ))
Local GridY:Int = Max(0,Floor( P.PositionY*TSPH.INV_SMOOTHING_LENGTH ))

seems to fix that.


Noobody(Posted 2009) [#21]
Thanks for the advice, DavidDC. I fixed this along with a few other minor bugs.

Also, I was able to find a new rendering method that is reasonably fast and still produces neat results:

Download (BMax source & Exe): Link

I used a method similar to the one described in Imaginary Human's article - still a prerendered Metaball texture, but a slightly different combination of blend functions.


In the past three weeks I also wrote on a raytracer in C++ (I began in BMax, but C++ is so much more comfortable due to operator overloading - vector/matrix calculations look much neater in C++ :) ).

The core of the raytracer is a fast triangle rendering method aimed at grid-based systems (127'000 trigs in about 7.2 seconds with a 800x600 resolution and 4x antialiasing). I then expanded the SPH code to work in three dimensions and wrote a small C interface for the BMax code to control the raytracer. The raytracer converts the SPH particles to a polygon using metaballs and the marching cubes algorithm. This gives very nice results, even though it's not realtime anymore.

I uploaded two videos on Youtube showing the SPH raytracer, here and here. I'm still working on this and am most probably going to upload a few more videos in the next weeks (with obstacles, more water etc.).


ImaginaryHuman(Posted 2009) [#22]
I like that it doesn't have any gaps in the larger areas of the fluid - and it is probably twice as fast as your other version. I peeked at the code, it looks pretty advanced and I couldn't understand it. :-)


Beaker(Posted 2009) [#23]
I have a problem running the code that comes with this (tho the exe runs fine). I get an error on this line:
ParticleTex = GLTexFromPixmap( LoadPixmapPNG( "Metaball.png" ) )
(actually the error is in this line:
If width=1 And height=1 RuntimeError "Unable to calculate tex size"
in GLAdjustTexSize which is called by GLTexFromPixmap)


The error is:
Unhandled Exception:Unable to calculate tex size

I can't really see why. Any ideas? Might be one for Mark to look at.

Hope this helps.


ker2x(Posted 2010) [#24]
You can speed up the Nice rendering by changing :
Const SAMPLE_RATE:Float From 3.0 to 5.0

(or more, but it's too ugly)

I'm trying to find a good candidate methode to do GPGPU using OpenCL.

While profiling the application to speed up the simulation part and rendering part, i found that a lot of time is not spent in Math but in memory transfert.

The slowest simulation Method is RedistributeGrid() which doesn't do much math, but a lot of memory transfert.

Before moving the math to OpenCL, the app need a major rewriting.


Hotshot2005(Posted 2010) [#25]
Impressive Stuffs :)


Taron(Posted 2010) [#26]
AH YES, I almost forgot about this, that was nice...
I took a look at it right now and replaced your metaball.png with a grayscale version from all the way white to black in the corners, then replaced your rendering loop with this:
			glBegin( GL_QUADS )
				Local S:Float = 8.0
				
				For Local P:TParticle = EachIn Particles
					Local col:Float = (0.5+0.1*P.Pressure)*P.Density
					glColor4f(col^3,col^1.5,col,1.0)
					glTexCoord2f( 0.0, 0.0 ); glVertex2f( P.ScreenX - S, P.ScreenY - S )
					glTexCoord2f( 1.0, 0.0 ); glVertex2f( P.ScreenX + S, P.ScreenY - S )
					glTexCoord2f( 1.0, 1.0 ); glVertex2f( P.ScreenX + S, P.ScreenY + S )
					glTexCoord2f( 0.0, 1.0 ); glVertex2f( P.ScreenX - S, P.ScreenY + S )
				Next
			glEnd()


That's kinda neat and no slowdown really.


ker2x(Posted 2010) [#27]
Change commited to http://github.com/ker2x/BM_Fluid
The standalone .exe can be found here : http://ker.endofinternet.net/img/XSPH.exe

Btw, you can find some of my old project in purebasic at http://ker.endofinternet.net/img/ :)

I'm planning to add GPGPU support to the simulation, using bah.opencl module : http://www.blitzbasic.com/Community/posts.php?topic=87452

I'm still waiting for a mail reply from the original coder to add the license in the code. For now, i assume it is released into some kind of opensource license.


ker2x(Posted 2010) [#28]
Mmm, i tought that the original coder was Tachyon, but it's noobody ?
Oops ... sending a new mail :)


smilertoo(Posted 2010) [#29]
Very nice effect, would be great if it could use gpu.


Taron(Posted 2010) [#30]
If I think about the troubles with openGL, I don't even wonna know what openCL is going to do...haha. But, yeah, wished it was straight forward and reliably compatible.


nawi(Posted 2010) [#31]
Beaker, try to change the texture to 2^n size.


ker2x(Posted 2010) [#32]
i find OpenCL easier than GLSL.
And you can do OpenCL without OpenGL.

You don't even need a graphic context to do OpenCL.
Or, better, you can do OpenCL without GPU :)


ker2x(Posted 2010) [#33]
I improved the rendering based on Metaballs.
Note the 556 FPS when the simulation is paused.
it run at 120FPS when unpaused.



The full source code can be found here : http://github.com/ker2x/BM_Fluid

I'm still rewriting the full app (and trying understand all the code).
I can't yet find a way to improve speed with OpenCL.


Nate the Great(Posted 2010) [#34]
I am quite jealous of your xsph skillz... more so of your rendering skills though. I am trying a marching squares algorithm and it is not going too well lol. I dont think the average gamer's computer is ready for something with as much action/particles as pixel junk shooter or even what I tried to do with flood but the frame rate got too choppy... sigh :/


ker2x(Posted 2010) [#35]
I'm not sure if you're replying to noobody or me. but... just in case :

Noobody is the original coder, i'm just (heavily) hacking it.
I still consider myself as very new to BlitzMax, never completed a full project with it. I'm not even a coder, i'm a sysadmin. :)

i don't understand (yet) most of the math behind it, noobody probably does but it's not even mandatory.

All you need is some passion and time, time, time, time (and github :p ) ... The more you're skilled, the less time you need. But nobody care/ask about how many hours/days/days/weeks/years were spent in this project. If you need 3 years of your spare time to write your marching squares algorithm ... well ... just do it. :)


Nate the Great(Posted 2010) [#36]
I was referring to noobody but thanks for the advice anyway. I replied like that because I wrote flood (link here) and that was my main problem.. rendering and graphics. I am redoing it now, hopefully much 'cleaner' this time.


ImaginaryHuman(Posted 2010) [#37]
I would think that probably the effort involved in marching cubes, or rather, marching squares for 2d, which btw is patented, is more than needed to simple draw some metaball images, especially if the metaball images are one or a few small images on a single texture, and especially if you are drawing using either vertex arrays or even point sprites. Vertex arrays at least should speed it up a lot. With metaball images you are just reliant mostly on fill rate. You should be able to do plenty of particles in realtime.


Nate the Great(Posted 2010) [#38]
@IH

yeah I have discovered that the hard way. and I read its patent expired in 2005. :) its not that its usefull for me now anyway...

and I honestly dont know how I am going to render metaballs... I would like to do it the same way Noobody has done it but I dont know what he would think and I honestly dont understand the rendering code.

edit: noobody, I raised the gravity to 2 and took off the 'speed limit' thing and the simulation explodes... with the speed limit on and high gravity it acts very odd still... just wondering if you know why this is?

edit2: still explodes with gravity 1 and speed limit on :/


ImaginaryHuman(Posted 2010) [#39]
See my tutorial in the tutorials forum about rendering metaballs using images. Once you've generated an `energy field` image you just draw it in LIGHTBLEND to add it to the scene. My guess is that there is some lightblending going on in the above image which causes the high-density/overlapping blobs to produce a whiter/lighter coloration.


Nate the Great(Posted 2010) [#40]
hey IH I saw your tutorial on it and that is how I did the graphics for my game flood. I just wish there was another faster less fuzzy looking way... but I am no good with gl ...


ker2x(Posted 2010) [#41]
I couldn't make it explode with the anti-cancer-spreading routine enabled.

If i comment the anti-cancer-cure it is, indeed, exploding at gravity 1.0


Noobody(Posted 2010) [#42]
If both the cancer cures are on (the density limiter and the speed constraint), it's impossible for cancer to spread. It is possible however, that the fluid behaves extremely unrealistic in an environment where normally cancer would blow up the fluid. This is due to the density limiter, which resets the density of cancer particles to a fixed value. This violates the conservation of mass. It's not that big of a problem if only one or two particles are affected, but if it's half of the fluid, the simulation looks really bad.

I really have no idea how I could fix this. All articles I read about this subject either prefer the gas pressure approach, which results in a compressible (read: unrealistic) fluid or use the approach proposed by Monaghan (which I am using) and implement some hacky routines to increase the stability.

The only idea I have right now is to run the simulation, remember particles that later become cancer and rerun the simulation, tracking those particles frame by frame to see what's going wrong.



Since I haven't posted in this thread for a while, I might as well talk about the progress that was made on the simulation. I'm working on a fluid module for BMax with support for multiple water tanks, different fluids with different properties, surface tension, soft bodies and rigid bodies.

I'm having a few problems with the rigid bodies, mainly because energy conversation is pretty difficult in a particle based approach :) But other than that, it comes along nicely. I made a small demo to demonstrate:

Screenshot:

Download (BMax code & exe): Link

Keys 1, 2 and 3 switch between the different examples. #3 is very, very WiP, since rigid bodies don't really work yet. In example 1 and 2, space triggers the fluid emitters. #1 uses the right and left mouse button for fluid/soft body placement and #3 only the left button for the rigid body creation.



I'm also working on my raytracer used to render the 3D version of the module. Currently I'm implementing photon mapping to allow the rendering of nice water caustics.

Color bleeding and caustics already work, but there are still a few issues, like the dark corners.

Here's an older video of the raytracer + 3D fluid in action.

So yeah, that's about it. I don't know if I have enough ambition to complete both the module and the raytracer, but I'll see where this is leading.



@NateTheGreat:
I am trying a marching squares algorithm and it is not going too well lol.

Marching squares on a screen-sized grid is, even for the GPU, an enormous overkill. My method is much simpler - it's basically the one in Imaginary Human's tutorial with a different blend mode:
You can increase the number of iterations in the loop at the end to reduce the fuzziness, but 2 iterations work pretty good. Full-screen blended quads take their toll on the performance, but that's the fastest method I could find. Another drawback is that only full colors can be used (R, G and B have to be either 255 or 0). This limits the amount of colors available to 8, which is... not good. I have no idea how to overcome that, so if anyone has suggestions - please tell me :)


ker2x(Posted 2010) [#43]
A good idea, i think, could be to transform overdense/overspeed/overheated water particle into vapor.

Nice raytracing... in BMax ? :)


ImaginaryHuman(Posted 2010) [#44]
You should at least implement vertex arrays, you'll probably see a 200% speedup.

Also maybe there are faster ways to calculate fluid, there seems to be a lot of work being done to make it happen. At least multithread it for a speedup? Also maybe switching to arrays instead of separate `TParticles` may be faster.


Nate the Great(Posted 2010) [#45]

If both the cancer cures are on (the density limiter and the speed constraint), it's impossible for cancer to spread. It is possible however, that the fluid behaves extremely unrealistic in an environment where normally cancer would blow up the fluid. This is due to the density limiter, which resets the density of cancer particles to a fixed value. This violates the conservation of mass. It's not that big of a problem if only one or two particles are affected, but if it's half of the fluid, the simulation looks really bad.

I really have no idea how I could fix this. All articles I read about this subject either prefer the gas pressure approach, which results in a compressible (read: unrealistic) fluid or use the approach proposed by Monaghan (which I am using) and implement some hacky routines to increase the stability.



I was nit picking lol ;) just tryin to see how much I could change and still keep it stable...

Marching squares on a screen-sized grid is, even for the GPU, an enormous overkill. My method is much simpler - it's basically the one in Imaginary Human's tutorial with a different blend mode:


yeah... just thought I would try it anyway.. and I am having a lot of trouble understanding opengl but thanks for the code

IH... I fail to see how implementing vertex arrays would be so much faster... seems like you would need so many of them it might be slower... ;) I dont know what im taking about


nice demo btw... Seems a bit buggy as far as water particles frequently getting stuck in the red stuff... but as you said, it is a wip. :)


ImaginaryHuman(Posted 2010) [#46]
Um... well in my tests using vertex arrays rather than immediate mode OpenGL primitives is usually about 2 times faster to render. There's significantly less function-call overhead. Also the use of arrays is faster than a linked list particle system.

Btw possibly a different curvature of your blobby image could allow the blob image to actually be smaller and give the same result? Ie a more intense curve might dropoff faster so that the image doesn't need to be so big - you could save on a lot of overdraw.

Another thought that came to mind is that if you're drawing fullscreen quads in order to `do math` on the rgb data, you might gain some speed if you break your screen into a grid of `dirty rects` and only draw small quads over the rectangles that actually have particles in them.


Noobody(Posted 2010) [#47]
A good idea, i think, could be to transform overdense/overspeed/overheated water particle into vapor.

Nice raytracing... in BMax ? :)

I once tried removing such particles (which is about the same as transforming them into vapor), but it actually makes things worse.

The raytracer is written in C++, because matrix/vector operations are so much more comfortable with operator overloading. Aside from that, the garbage collector in BMax would have taken a huge toll on the performance due to the amounts of (vector) objects used. I once wrote a test implementation of a raytracer in BMax but dropped it real quick in favor of C++ because the speed was just inacceptable. I still use BMax to calculate the fluid and to control the raytracer, though.

You should at least implement vertex arrays, you'll probably see a 200% speedup.

Also maybe there are faster ways to calculate fluid, there seems to be a lot of work being done to make it happen. At least multithread it for a speedup? Also maybe switching to arrays instead of separate `TParticles` may be faster.

Vertex buffers would be a good idea, the rendering method needs a makeover anyway :)

The fluid calculation is pretty much optimized. Values are precalculated where possible, a grid/linked list combination is used to cull unimportant particles and divisions are avoided as much as possible. Multithreading is somewhat planned, but arrays instead of types are not an option. It may be faster, but the code would look a lot uglier and more confusing.


ker2x(Posted 2010) [#48]
Multithreading is somewhat planned, but arrays instead of types are not an option. It may be faster, but the code would look a lot uglier and more confusing.


I will have to do that if i want to write the openCL implementation of XSPH. But it's not planned for the next few days.


ZJP(Posted 2010) [#49]
Hi,

W.I.P for Blitz3D and Unity http://www.zinfo972.net/avatar/FluidDEMOV4.zip

]


JP


Noobody(Posted 2010) [#50]
3D update!




ImaginaryHuman(Posted 2010) [#51]
Hey very nice, man. It doesn't quite look like water to me, maybe it's the particle size, looks a bit more like a gloopy thicker solution, but very nice rendering.


Noobody(Posted 2010) [#52]
Thanks!

Yeah, I'm still trying to figure out why it looks that way. I suspect it's because it moves too slow and therefore looks as if it had a high viscosity.

I'm going to try skipping every second or third frame to make the video faster.


Hotshot2005(Posted 2010) [#53]
it is amazing on what you have done! :)


*(Posted 2010) [#54]
another thing is the stuff coming out the pipe looks like glue there is no activity on the water coming from it where normal water does something.


Noobody(Posted 2010) [#55]
Yeah, that's because the particles in the emitter are created in a regular grid, which is why the water-jet always keeps the same shape. Adding a random offset to the particles when created at the emitter should fix that.


Kryzon(Posted 2010) [#56]
Very nice refractive material!

but the fluidness needs some work, and I don't think that's related to the video speed (even slow motion water should look watery).

I'm thinking decreasing the "weight" of each particle should improve it. Also, if there is any friction between each particle, you should reeaally decrease it: takes too much for the bounced wave to get back to the flow of water. You can see from 00:05 up to 00:16 that the water takes too much time to cover itself, and this sort of "depression" appears (even if the impact of the water-jet has some influence on that, it still doesn't feel natural).


Rck(Posted 2010) [#57]
great to see the fluid and raytracer working well together, makes me want to add photon mapping to my own.


Noobody(Posted 2010) [#58]
Today I experimented with multithreading to speed up the SPH algorithm. Unfortunately, the changed GC behaviour in threaded build is much slower, which results in more speed lost than could be gained through multiple threads (when using many objects, that is).

For comparison, I ported the code to C++ and implemented multithreading directly using the WinAPI. Not only did the program (without threading) gain 20% compared to the BMax version, but multiple threads actually increased the speed.

This means that the best way probably is to do most of the work in C++ and only do rendering and initialisation in BMax. Getting threading to work in C++ on all three platforms requires additional platform specific code and a lot of testing, but the speed is worth it, I guess :)

I compiled a small realtime 3D demo of the C++ SPH project to showcase the new multithreading option:

Download: Link

The cube can be turned using the arrow keys. Since the gravity always pulls downwards, it's a good way to play around with the fluid. On startup, the program asks how many threads it should use - just input the number of cores in your CPU for optimal performance, but not more than 4. With 4 threads, the performance peaks. The FPS gained trough more than 4 threads are lost due to bigger overhead.

Next I'll be 'multithread-ising' the raytracing code to decrease the rendering time so I don't have to wait hours for a video to finish :)
After that I'll probably work on a SPH module implementing the 2D multithreading version in C++.


but the fluidness needs some work, and I don't think that's related to the video speed (even slow motion water should look watery).

I *think* it's because of the viscosity constants that were set to high. I fixed it in the demo above and I think it's better now. I'll try to render a new video to see whether it helped. Thanks for the advice!

great to see the fluid and raytracer working well together, makes me want to add photon mapping to my own.

For implementing photon mapping I can recommend Jensen's book "Realistic Image Synthesis Using Photon Mapping". It explains the algorithm in great detail and provides a demo implementation in C++. The math he uses is mind-bending, but when it comes down to code, it's pretty straightforward. For diffuse interreflection, photon mapping is too noisy in my opinion, but it yields great results for caustics.


slenkar(Posted 2010) [#59]
i wonder how much faster it would be in plain C?
:)


Robert Cummings(Posted 2010) [#60]
Slower would be my guess? :P


Noobody(Posted 2010) [#61]
New video!



More photons, more particles, more threads, more everything! Except render time, of course.

Next stop: obstacles.


Clyde(Posted 2014) [#62]
Well Done!! :)

Very very cool indeed!!